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Chapter 1 A) 
Introduction to GeomInt2 zur: 


Thomas Nagel, Tuanny Cajuhi, Ralf-Michael Günther, Jobst Maßmann, 
Frank Wuttke, Keita Yoshioka, and Olaf Kolditz 


1.1 Project 


GeomInt2 is the follow-up project to the joint project “Geomechanical Integrity of 
Host and Barrier Rocks—Experiment, Modelling and Analysis of Discontinuities 
(GeomInt)”. The results of GeomInt made essential contributions to the analysis 
of the origin and development of discontinuities in clay, salt and crystalline rocks. 
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Fig. 1.1 General project 
structure 


wP 3 
WP 4 (Digitalization) (Crystalline) 


Discontinuities are considered in many forms, be it as volumetrically distributed 
damage macroscopically representing fissures in the rock microstructure, disconti- 
nuities that can form anew at phase interfaces as well as discrete cracks or fracture 
networks. The pathways created or extended by these discontinuities entail the risk 
of migration of fluid phases from deep to near-surface geological layers and must 
therefore be taken into account in geotechnical safety analyses. 

The main result of the primarily methodology-oriented joint project GeomInt is the 
development of a broad experimental and numerical platform for laboratory and field 
analysis, and prediction of discontinuities in important reservoir and barrier rocks. 
Model comparisons for damage and fracture processes driven by different physical 
phenomena provide indications for the optimal application areas of the numerical 
methods. The follow-up project GeomInt2 aims to demonstrate the practical appli- 
cability of the developed experimental-numerical instruments under real conditions 
at the site scale. In addition, knowledge gaps that still exist will be closed through in- 
depth tests on a smaller scale and in-depth model developments for specific process 
components. 

At the core of GeomInt2 was an interdisciplinary consortium of partners from uni- 
versities, state and private research institutions with complementary, long-standing 
experience in the analysis of geosystems that has proven itself in the previous project. 
New insights were elaborated, especially in understanding the effects of discontinu- 
ities on underground geosystems. In this context, the focus of experimental investiga- 
tions shifted significantly to the implementation and evaluation of in-situ experiments 
in various underground research laboratories (URLs). For the efficient numerical sim- 
ulation of coupled mechanical, thermal and hydraulic processes in the formation and 
development of discontinuities on the scale under consideration, the models and algo- 
rithms were further refined in GeomInt2. To enable the analysis of larger and more 
complex systems, high-performance computing (HPC) methods were added as a new 
aspect. The descriptive presentation of structural, experimental and model results in 
a real context (e.g. URLs) is to take place within the framework of an integrated 
visual data analysis approach (virtualisation). Figure 1.1 is a graphic illustration of 
the GeomInt2 project. 

The project results aim at improved understanding of the processes, the used meth- 
ods and the application-oriented systems for realistic time and length scales in order 
to make the planning and realisation of geotechnical uses of the subsurface safer, 
more reliable and more efficient. Another important advantage of GeomInt2 is the 
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Fig. 1.2 GeomInt book 


GeomInt—Mechanical 
Integrity of 


Host Rocks 


transferability of the experimental-numerical concepts and methods to other geotech- 
nical applications such as deep geothermal energy, energy storage, repository prob- 
lems, methods for hydraulic stimulation, conventional and unconventional resource 
extraction or tunnelling. GeomInt2 intensified the internationalisation started in the 
previous project by cooperation with complementary research projects (e.g. DECO- 
VALEX 2023 and the Mont Terri project) (Fig. 1.2). 

The work packages mentioned in Fig. 1.1 were closely linked to the URLs Mont 
Terri (WP1/WP2), Springen (WP2) and Reiche Zeche (WP3) in continuation of the 
first project phase. The sample material for the geotechnical laboratories in Freiberg 
(TUBAF), Kiel (CAU) and Leipzig (IfG) originated from these URLs. Data from in- 
situ experiments in Mont Terri (BGR) and the Springen pit (IfG) were incorporated 
into the modelling platform of the GeomInt project coordinated by the UFZ, in which 
all project partners participated. So-called Model EXperiments (MEX) were defined 
as a synthesis tool, in which experimental and modelling work were systematically 
combined from the beginning of phase one (Kolditz et al. 202 1a, b). This integration 
formed the basis for the continuation of the project and the application of its results 
to site-scale problems in phase 2. 


1.2 Approach 


The research approach of GeomInt2 continues following the general concept of 
discontinuities in geosystems as illustrated in Fig. 1.3. Discontinuities in geological 
media—being a general feature in various host rocks such as salt, clay and crystalline 
rocks—occur at various scales. New fractures can be initiated by elevated fluid pres- 
sures, temperatures, and/or mechanical stresses in conjunction with geotechnical 
applications (TH?MCB drivers Fig. 1.3). Existing joints or faults can be reactivated 
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Fig. 1.3. Graphical abstract of the GeomInt project: geological reservoir-barrier systems and geome- 
chanical integrity (Nagel et al. 2021) 


by similar multi-physical reasons as well (TH*MCB evolution Fig. 1.3). Fractures 
are localized and subject to size effects—they thus render a problem inherently 
multi-scale. Changes in connectivity of existing fracture networks due to evolving 
discontinuities can qualitatively alter system behaviour and lead to accelerated path- 
ways from repositories to cap or side rocks (near-to-far-field bridging in Fig. 1.3). 

GeomInt2 leveraged both the experimental and numerical framework from the 
Lab to Field Scale—which became the subtitle of the present book. GeomInt2 was 
aimed at improving our understanding of the underlying physical mechanisms and 
provide methods across different geological settings. The numerical instrumentation 
has been further developed and tested against experimental results on the field scale 
in particular from the underground research labs and mines in Mont Terri (CD-A 
experiment, Opalinus clay), Reiche Zeche (crystalline rock, STIMTEC experiment), 
and Springen (salt rock, large borehole experiment). 

As a result of the GeomInt project so far, a broad experimental and numerical 
platform has been developed for the investigation of discontinuities due to swelling 
and shrinking processes (WP1), pressure-driven percolation (WP2) and stress redis- 
tribution (WP3). With these phenomena, processes relevant in the most important 
barrier and reservoir rocks (clay, salt, crystalline) are addressed. A comprehensive 
validation of the platforms (,,Proof-of-Concept*) was done by the Model EXperi- 
ments (MEX) for the damage and fracture processes driven by different processes. In 
the continuation of the project, the basic idea—the investigation of the three relevant 
process types for the formation of pathways in different rock types—is preserved. 
The proven project structure (WP1-3) is also retained (Fig. 1.4a). 
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(a) Networking within work packages (b) Networking with STIMTEC 


Fig. 1.4 GeomInt project geography 


The emphasis in GeomInt2 is set as follows: 


— Experiments: Focus on the evaluation of recent in situ experiments at the Mont Terri 
(Opalinus clay), Springen (salt) and Reiche Zeche (crystalline) rock laboratories. 
Selected laboratory experiments are used to specifically fill remaining knowledge 
gaps (e.g. mechanical anisotropy of clay rocks). 

— Modelling: Numerical methods are to be further developed and used for in-situ 
applications. For the scaling from micro- to macro-scale models, different model 
approaches are to be specifically linked with each other (e.g. integration of the 
LEM into the FE method). 

— Digitalisation: New developments in the computer sciences, such as 
high-performance computing and visual data analysis, are to be incorporated into 
the continuation of the project. This is necessary, on the one hand, to achieve the 
required scaling of the models to the application scale and, on the other hand, to 
allow a presentation of the model results in the real context of the underground 
laboratories. 

— Internationalisation: In addition, the international cooperation with the Mont Terri 
project and DECOVALEX 2023 should be intensified, also in order to strengthen 
the external impact of the Geo:N project. 


Figure 1.4 depicts the project geography of GeomInt2 in direct continuation of 
the workflows established in phase one of the project between the project partners 
from the universities of Freiberg (TUBAF), Kiel (CAU) and Stuttgart (UoS) as well 
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as research institutions BGR, UFZ and the IfG company. The close collaboration 
established between the GeomInt and STIMTEC projects both working on disconti- 
nuities in crystalline rocks linked to the URL “Reiche Zeche” was also continued in 
the second phase. 

The following chapters introduce some important results from the three work 
packages at the core of the project: WP1 (swelling/shrinkage), WP2 (fluid percola- 
tion) and WP3 (stress redistribution) for three types of host rocks, namely, clay, salt 
and crystalline rocks. Subsequently an introduction to the virtual reality and compu- 
tation platforms is given, featuring the visualization, HPC and software engineering 
aspects. Some results are given more emphasis than others. Particularly, simulations 
on the hydro-mechanical behaviour of Opalinus clay were presented in greater detail. 
In the other chapters, a shorter presentation was chosen to keep the book concise. 
However, we point the reader to publications where more information also on the 
remaining topics addressed in the project can be found. The book closes with a 
synthesis and outlook. 
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Chapter 2 R) 
Hydro-Mechanical Effects and Cracking u 
in Opalinus Clay 


Tuanny Cajuhi®, Nima Haghighat, Jobst Maßmann, Mostafa Mollaali, 
Amir S. Sattari, Vahid Ziaei-Rad, Gesa Ziefle, Thomas Nagel, Frank Wuttke, 
and Keita Yoshioka 


2.1 Objectives and Outline 


In this chapter, we investigate hydro-mechanical effects in the Opalinus Clay, espe- 
cially those leading to cracking. We present a methodology comprising laboratory 
and field scale experiments, as well as the development and application of numerical 
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approaches. In Sect. 2.2, we present a short overview concerning the Opalinus Clay 
(OPA) formation, its behavior and applications. This potential host rock for radioac- 
tive waste disposal shows complex coupled behavior, especially related to shrinkage 
and shrinkage/desiccation-induced cracking. 

Section 2.3 describes the Cyclic-Deformation (CD-A) experiment in the Mont 
Terri Rock laboratory, one of the focal points of the GeomInt2 project. The experiment 
allows the measurement of relative humidity, deformation and cracking at in-situ 
scale. Furthermore, a drilling campaign allowed to extract core material that has 
been used in different laboratory experiments described in this chapter. We monitor 
the fracture response of OPA and evaluate the influence of its anisotropy (bedding 
orientation) as well as the effect of different suctions on the load-displacement curves. 
These experiments are modeled within numerical approaches based on the methods 
described in Sect. 2.4. In the latter, a continuum based, unsaturated hydro-mechanical 
approach is extended with the phase-field method (Sect. 2.4.2) and with the discrete 
element method (Sect. 2.4.3). 

After comparing the response of the material both at the laboratory and at the 
field scales, we propose an interactive visualization strategy in Sect. 2.6. Finally, 
we review the proposed methodology in Sect. 2.7 and discuss its extension, future 
applications and transfer. 


2.2 Observations and Processes in Opalinus Clay 


The Opalinus Clay formation (OPA) is considered as a potential host rock for radioac- 
tive waste disposal due to its low hydraulic permeability and ability to retain radionu- 
clides. Opalinus Clay formations are located in Northern Switzerland (Nagra 2021) 
and South-Eastern Germany (Hoth et al. 2007; Bundesgesellschaft fiir Endlagerung 
(BGE) 2020). In comparison to other clay formations, the OPA is well characterized 
by numerous in-situ experiments, conducted in the framework of the international 
Mont Terri Rock Laboratory in Switzerland. Its applicability in various geotechnical 
applications including the disposal of radioactive waste have been investigated, e.g. 
Bossart et al. (2018) and references therein. 

The formation was accumulated during the dogger period (Jura) sandy and sandy 
facies as well as related subfacies types (Lauper et al. 2018; Kneuker et al. 2021). OPA 
is an anisotropic and highly heterogeneous material and shows complex and cou- 
pled thermo-hydromechanical (THM) behavior, including shrinkage and swelling. 
Swelling is related to an expansion of the material. In the current study, we focus 
on shrinkage, which is related to bulk volume reduction induced by environmental 
changes such as variation of relative humidity (RH), that, in turn is related to changes 
in water content. This effect can be observed and measured both in laboratory and 
field scales (Wild et al. 2017; Bock 2000; Zhang et al. 2010). Furthermore, it can lead 
to shrinkage-induced or desiccation cracks. With focus on how such cracks develop, 
the behavior of OPA is investigated in the following via experimental and numeri- 
cal approaches that are applied at different scales. In-situ monitoring and borehole 
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extraction are carried out in the Mont Terri Rock Laboratory in order to characterize 
the THM-behaviour of the OPA (Bossart et al. 2018; Kneuker and Furche 2021). In 
particular, the CD (Cyclic Deformation) experiment was a valuable starting point to 
gather experimental data and experience at the in-situ scale and to investigate the 
hydraulic-mechanical coupling induced by swelling and shrinking of shaly facies of 
the OPA due to cyclic variations of air humidity (Matray and Mori 2012; Ziefle et al. 
2017). In the follow-up project, the CD-A-experiment, these processes were inves- 
tigated further in detail, this time in the sandy facies of the OPA. Consequently, the 
CD-A-experiment is one focal point of the GeomInt2 project and will be described 
in the following. 


2.3 Monitoring and Data Integration 


In the Cyclic-Deformation (CD-A) experiment (Ziefle et al. 2022c), coupled effects 
due to the excavation and seasonal changes in relative humidity (RH) are monitored 
regularly. The long-term measuring program focuses on a detailed characterization of 
the host rock as well as on the evolution of process variables like the pore pressure, the 
displacements and the water content amongst others. The evolution of the excavation 
damaged zone (EDZ), corresponding changes of the water content and monitoring of 
desiccation cracks within the first decimeters of the rock are the central points of the 
experiment. These measurements enable an increased understanding of coupled HM 
effects in the EDZ and provide a valuable basis for a multi-disciplinary investigation 
aiming at further development of numerical modelling approaches and constitutive 
equations to reliability capture coupled effects in the OPA. 

The CD-A experiment is characterized by two horizontal niches (twins) that are 
excavated perpendicular to the strike of bedding and that are not covered with any 
kind of stabilization (like shotcrete). Consequently, a direct interaction between the 
climatic conditions and the host rock takes place. While the HM effects in the open 
twin are evoked by the excavation and the seasonal variation of temperature and 
relative humidity, the relative humidity in the closed twin remains constant due to 
a sealing at its entrance. With that, the desaturation effect is minimized as much as 
possible holding the relative humidity high. The different evolution of the measured 
relative humidity within the niches is pictured in Fig. 2.1. The sudden drops in the red 
curve (closed niche) are related to door openings and installations. The open niche 
shows shows a cyclic decrease and increase of RH correspondent to the winter and 
summer months, respectively. 

The changes in RH, especially those leading to an increase of suction, i.e., decrease 
in the capillary pressure, can lead to cracking. This phenomenon has been observed 
in both twins, but more intensely in the open niche. A detailed monitoring of the 
evolution of desiccation cracks along a horizontal line (scanline) in each niche is 
presented in Regard et al. (2021). The program includes monthly monitoring of the 
length and the aperture of the cracks initiated at the niche walls. The measurements 
have indicated significantly more desiccation cracks in the open than in the closed 
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Fig. 2.1 CD-A experiment: measured relative humidity [%] in the twin niches. For example, the 
period between September 2019 and February 2020 shows a decrease of the RH, i.e. desaturation. 
Note that the RH in the closed niche remains high and approximately constant 


twin. The sum of crack apertures in the open niche is also larger than in its closed 
twin. Due to the increase of relative humidity in the summer months, the sum of the 
crack apertures decreases. 

In addition to the RH and cracking measurements, the CD-A experimental cam- 
paign comprises measurements focusing on the characterization of the claystone 
as well as the observation of process variables like the convergence, pore pressure 
and permeability which are amongst others measured in the boreholes presented in 
Fig. 2.2. Furthermore, laboratory scale tests are part of the study and provide further 
understanding on the OPA behavior. These tests have been carried out at material 
from CD-A drill cores extracted in March 2021 in the context of the permeabil- 
ity measuring campaign as presented in Fig. 2.2. Multi-disciplinary investigations 
along the presented boreholes are discussed in Ziefle et al. (2022a) and underline 
the high heterogeneity of the host rock with different facies and sub-facies types and 
the impact of the excavation process and the ventilation conditions on the near field 
around the niches. 

In the following, we describe the performed laboratory experiments that focus on 
anisotropy characterization and mechanical response of the material under different 
suctions. 


2.3.1 Laboratory Experiments 


2.3.1.1 Mechanical Characterization of Opalinus Clay 


In the first phase of the GeomInt project, fracture propagation characteristics in 
Opalinus clay have been investigated both numerically and experimentally. In this 
regard, a set of three-point-bending experiments has been performed on rectangular 
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Fig. 2.2 Summary of boreholes from measuring campaign in 2022. Extensometers, mini- 
piezometers and permeability sensors are installed, among others 


samples of Opalinus clay. Given the status of bedding structure in obtained cores, 
only the configuration of loading being perpendicular to the bedding direction could 
be studied. Accordingly, the effect of structural anisotropy on sample strength as well 
as the fracture propagation behavior was investigated, albeit qualitatively, adopting 
different numerical methods. As an extension to the first phase of the project, a series 
of semi-circular bending experiments has been conducted in order to thoroughly 
investigate the direction-dependent strength characteristics of the Opalinus clay. The 
semi-circular configuration (Fig. 2.3) was chosen because of its minimal material 
requirements and relatively easier preparation procedure. The samples were taken 
from two different drill cores of sandy facies type, namely A24 and A26, near the 
CD-A twins (Kneuker and Furche 2021; Ziefle et al. 2022a,c). 

Following ISRM suggested regulations (Kuruppu et al. 2014), the samples were 
prepared with radius of 36 x 107° m, thickness of 40 x 107? mand an artificial notch 
at center with length of 14 x 107° m. Given sensitivity of the clay material to water, 
all specimens were prepared under dry conditions using a diamond band saw with 
2 x 10-7 m thickness. As summarized in Fig. 2.3, three type of samples, namely P-, 
Z- and S-samples, have been considered for the experiments. The naming convention 
is indicative of loading angle 0, which represents the orientation between the isotropy 
plane and the loading direction. We have tested following loading directions 6 = 
[0, 45, 90]° that represent P-, Z- and S-samples, respectively. 

In total, three samples, two from A26 cores and one from A24 core, for each 
loading angle were tested. The tests were performed under quasi-static loading con- 
ditions with strain rate of 0.03 mm/min. The typical laboratory results in terms of 
applied load versus the deflection of the samples obtained from A24 core drill are 
illustrated in Fig. 2.4a. The load-deflection curves in the three cases exhibited a sim- 
ilar pattern; a downward convex trend followed by an approximately linear increase 
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Fig. 2.3 Schematic summary of SCB samples: (from left to right) P-samples, Z-samples and 
S-samples 


Load [KN] 


Peak Load [KN] 
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Fig. 2.4 a Obtained Load-Deflection curve for Opalinus clay sample with different loading angles, 
b Mean value and range of the obtained peak loads for each loading angle. With increasing loading 
angle, an increase of the peak load is observed 


in the load response until failure. Furthermore, the sharp drop of the load repre- 
sents the brittle failure mechanism during the experiments. The mean value and the 
range between the maximum and minimum achieved failure loads are summarized in 
Fig. 2.4b. As shown, the maximum obtained load before failure, indicative of mate- 
rial toughness, increases as the angle between the loading direction and the isotropy 
plane increases. The considerable range among toughness of the samples with sim- 
ilar anisotropy structure could be attributed to various factors, e.g. micro defects 
generated in samples during the preparation process or the conditions at which cores 
were drilled. 

Figure 2.5 shows the final fracture pattern of the three samples from A24 cores. 
Note that the crack does not propagate straight along the bedding direction for the 
P-sample due to the heterogeneity of the material. In the Z-sample, the crack prop- 
agation along bedding direction can be identified more clearly. Finally, the crack 
propagation in the S-sample also shows, initially, a strong dependency on the bed- 
ding direction, but continues its trajectory straight, as expected. 
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Fig. 2.5 Failure pattern and fracture path trajectory of a P-sample, b Z-sample and ce S-sample. 
The cracks propagate approximately straight from the notch 


2.3.1.2 OPA Characterization Under Different Humidity Conditions 


As also observed at in-situ scale, in the CD-A experiment, changes in the RH lead to 
variations in the mechanical and fracture behaviors. We investigate this phenomenon 
in the laboratory by characterizing the mechanical response of the OPA at various 
RH. To do that, a customized desiccator cell has been configured and set up at the 
geomechanics and geotechnics laboratory of CAU. The idea originates from the first 
phase of the project, where the concept of a humidity controlled bending test setup 
was proposed (Sattari et al. 2021). As shown in Fig. 2.6, the customized desiccator 
cell contains four distinct parts, namely 


1. Epoxy glass cell, equipped with RH and temperature sensors, to provide an 
isolated environment for the designed experiments 

2. Salt bath 

3. Interchangeable mechanical loading apparatus and 

4. Adjustable camera tube fabricated for use of a micro-camera dedicated to optical 
displacement measurements. 


The test procedure generally takes place in two steps; the equilibrium phase and the 
loading step. The sample is exposed to a certain RH during the equilibrium phase 
while being allowed to adjust with the environment. The humidity conditions inside 
the cell were achieved and controlled by applying different supersaturated salt solu- 
tions. For the Opalinus clay samples, the equilibrium normally occurs between 7 to 
10 days. The equilibrium period varies due to different initial water content of the 
sample and the specified RH of the cell. The equilibrium conditions are monitored 
with consecutive weight measurements of each sample. The unique design of the 
cell provides an opportunity to assess mechanical characteristics of OPA under con- 
trolled humidity conditions with minimizing the environmental exposure of sample 
between equilibrium phase and loading step. In this section, anisotropic strength 
characteristics of Opalinus clay are studied in different humidity conditions using 
the semi-circular bending configuration. Table 2.1 summarizes the four different 
humidity conditions considered in this study. The obtained suction pressure for each 
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Fig. 2.6 Customized Desiccator Cell equipped with relative humidity (RH) and temperature sen- 
sors, a salt bath, an interchangeable mechanical loading apparatus and an adjustable camera tube 


Table 2.1 Utilized salt solutions along with the achieved relative humidity and generated suction 


Salt solution Temperature (°C) Relative humidity Suction (MPa) 
KNO3 23.4-24.5 0.88 17.5-17.55 
KCl 21.2-22.5 0.82 26.6-27.07 
NaCl 21.2-22.5 0.75 38.4-39.6 
Room conditions 22.5-23.0 0.42 118.3-118.5 
LiCl 24.5-25 0.19 217.6-229.6 


salt solution using 3 samples each is estimated using Kelvin’s equation and based on 
the recorded RH and temperature during the experiments, 


pee i (2) (2.1) 
Vino M Po l 


where the universal gas constant R = 8.314 J/(mol K), T is the temperature in Kelvin, 
Vwo is the specific volume of water and M = 0.01801528 kg/mol is the molar mass 
of the pore fluid. The term p/ po corresponds to the relative humidity (RH). 

The obtained experimental results are presented in terms of failure load as a func- 
tion of imposed suction for each loading angle in Fig. 2.7. A general increasing trend 
for the material strength following the increase in the imposed suction and resultant 
decrease of water content can be clearly observed. This behavior is associated with 
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Fig. 2.7 Load-displacement curves under different RH. In general, the load increases with suction 


the generation of capillary forces in a partially saturated porous medium due to the 
pressure difference between gas and fluid phases at the pores (Wild et al. 2015). 

Considering the aspects mentioned in this section, we aim at reproducing the 
shrinkage-behavior of OPA numerically in the following. 


2.3.1.3 Additional Tests on Percolation 


We utilized the drilled cores to carry out fluid percolation tests. These tests are 
important in the context of hydraulic fracturing. Knowing when the fluid will induce 
a crack aids in understanding critical conditions in the field. Fluid percolation exper- 
iments were carried out in the geomechanics and geotechnics laboratory of CAU 
Kiel. The main objective to these experiments is to investigate the characteristics of 
fluid driven fractures in Opalinus clay under effects of structural anisotropy as well 
as the anisotropy of the stress state. The utilized setup includes a true triaxial loading 
apparatus, resembling the in-situ stress conditions, and a syringe pump. More details 
related to experimental setup and the true-triaxial device can be found in Sattari et al. 
(2021). Cubic samples were prepared with length of 43 x 107°? m and a scaled bore- 
hole drilled with diameter of 7 x 107° m and length of 20 x 107° m. Fluid is injected 
under constant flow rate into the sample through an adapter and a steel tube glued 
to the drilled borehole. A fluorescent oil was used as the injection fluid in order to 
trace even the tiniest fractures initiated during the process. Schematic representation 
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Fig. 2.3 Schematic representation of cubic S-, P- and Z-samples (left to right) 


Table 2.2 Summary of fluid percolation experiments on Opalinus clay 


Stress configuration (MPa) P sample S sample Z sample 
Sy = 12, Sy = 12, Sy = 12 v v v 
Sy = 12, Sy = 8, Sh = 8 v vV - 
Sy = 12, Sy = 8, Sh = 4 vV v = 


of the applied principle stresses along with the bedding direction of the samples are 
shown in Fig. 2.8. 

Given the position of the scaled borehole with respect to the bedding direction, 
three different loading configurations can be tested: 


1. S-samples where the drilled borehole is perpendicular to the bedding direction 
2. P-samples where the borehole is parallel to bedding and 
3. Z-samples where bedding direction is inclined to the borehole. 


Accordingly, Table 2.2 summarizes the number of tested samples for each bedding 
direction and the stress state in which the experiment was conducted. 

The experimental test starts after the in-situ stresses are applied. Once the mechan- 
ical equilibrium has been reached, the fluid injection then starts at flow rate of 10 ml/h. 
Upon observing the breakdown pressure, the experiment is stopped and the sample 
is retrieved. After careful investigation of failure patterns on the surface, the samples 
were sawed into number of cross section to have a better understanding of frac- 
ture propagation behavior inside it. Figure 2.9 shows the crack morphology on the 
surface, the outline of the cross sections, and the inner trajectories of fracture path 
indicated by florescent oil under UV light. It can be clearly seen that the pressure 
induced fractures inside Opalinus clay can either follow a straight path or propagate 
in branches and form a network. These observations are also reflected in the pressure 
response shown in Fig. 2.10. Generally, the drop in the pressure response is associ- 
ated with increase in the occupied volume by the injection fluid followed by initiation 
and propagation of fractures inside the sample. In some cases, e.g., Fig. 2.10c, the 
initiated fracture reaches the surface of the sample rapidly and causes the depletion 
of the injection fluid out of the sample and a sharp breakdown in pressure response. 
In other cases, e.g., Fig. 2.10b, fractures propagate in tiny branches and form a net- 
work before reaching the sample surface which is reflected in the fluctuations of the 
pressure response. 
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Fig. 2.9 Fracture path trajectory, cross-section cut and inner fracture path trajectory indicated by 
florescent oil for a P-sample b S-sample and e Z-sample of Opalinus clay under isotropic in-situ 
stress of 12MPa 
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Fig. 2.10 Comparison between pressure-time curve for different isotropy planes of Opalinus clay 
under mechanical stress configurations of a Sy = 12 MPa, Sy = 12 MPa, Sh = 12 MPa, b Sy = 12 
MPa, Sy = 8 MPa, Sh = 8 MPa, and c Sy = 12 MPa, Sy = 8 MPa, Sh = 4 MPa 
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2.4 Numerical Modeling Approach 


To reproduce shrinkage effects observed in OPA in the laboratory and field scales 
(CD-A experiment), we utilize and extend the numerical methods evaluated in the 
first phase of the GeomInt project (Yoshioka et al. 2021; Kolditz et al. 2021a). The 
extensions are related to the anisotropic and cracking behaviors of the material. 

Intuitively, cracks generate edge separations with creation of boundary surfaces 
over which the displacements are discontinuous. This makes numerical simulation of 
crack propagation challenging especially if we do not prescribe the location, topol- 
ogy, or time of their creation. For modeling crack propagation, several strategies are 
proposed in the literature. One way to deal with this sort of topology changes is by 
adaptively remeshing the domain (Bouchard et al. 2000; Meyer et al. 2004; Branco 
et al. 2015; Rangarajan et al. 2015). However, assigning boundary conditions to the 
newly created surfaces and including potential contact are not straightforward tasks. 
Another way to model crack propagation is to enrich elements that contain cracks to 
capture the displacement discontinuity, an approach known as eXtended/Generalized 
Finite Element Methods (XFEM, GFEM) (Moés et al. 1999; Belytschko et al. 2001; 
Gasser and Holzapfel 2005; Fries and Belytschko 2006; Belytschko et al. 2009; 
Khoei et al. 2012). Despite recent success (Shauer and Duarte 2022), the implemen- 
tation efforts remain a daunting task especially when arbitrary fracture topologies 
are considered (Belytschko et al. 2001, 2003; Duarte et al. 2007). 

Crack initiation needs to be prescribed in the methods discussed above. This can 
become a limitation in applications of geomechanical integrity, where we need to 
assess fracture initiation possibilities without any a-priori knowledge. The phase- 
field approach for brittle fracture represents the crack using a phase-field (damage) 
variable that results in a diffuse cracked interface. Furthermore, it does not require any 
special mesh treatment and prescription of the crack position. We apply and evaluate 
the response of both diffuse cracks using a continuum mechanics based approach 
combined with the phase-field and sharp cracks using the discrete element method. 
While the former is based on successive minimization of the total energy and cracks 
initiate as a result of the energy minimization, the latter finds fracture initiation based 
on the stress as the stress exceeds the strength defined between discrete elements. 


2.4.1 Unsaturated Hydro-Mechanical Approach 


The continuum formulation of the unsaturated HM-coupled problem outlined below 
is discussed in detail in Pitz et al. (2022) and is mainly based on a simplification of 
the approach from Grunwald et al. (2022). It is implemented in the current version 
of the open source code OpenGeoSys (Bilke et al. 2022). 

The isothermal formulation is founded on the conservation of water volume in 
a deformable porous medium according to Biot’s consolidation theory (Biot 1941) 
considering unsaturated Darcy flow by the Richards equation (Richards 1931): 
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with the density of the liquid phase prr, the liquid saturation SL, the porosity &, the 
liquid and solid compressibility £p Lr, B),sr, respectively, the Biot coefficient ag, 
the liquid pressure pyr, the intrinsic permeability tensor ks and relative permeability 
ke the liquid viscosity urr and the solid displacement vector u. 

The capillary pressure p..p is defined as the difference between the liquid pressure 
pır and the gas pressure por. Following the Richards equation, pressure gradients 
in the gas phase are neglected and by choosing the gas pressure to zero, we obtain 
following relationship 


Pcap = PGR — PLR = —PLR- (2.3) 
The liquid saturation SL is linked to the effective liquid saturation via SL = 


Se (St sna _ Shires) + SL,res, Which is in turn defined as a function of the capillary 
pressure according to van Genuchten (Van Genuchten 1980): 


n 1-1 
S.= (1 ve (=) ) , (2.4) 
Pe 


with pe as the gas entry pressure and n € [0, 1) as empirical material property. The 
relative permeability is defined according to van Genuchten/Mualem (Mualem 1976): 


=, o> 


The momentum balance equation is given by 


0 = div (o’ — æg SL piel) + [(1 — $) Psr + Sibpirlg. (2.6) 


where psr is the density of the solid phase, g the earth’s gravitational accelera- 
tion vector and o’ represents the effective stress tensor, which can be calculated by 
the solid stiffness C and the elastic strains: o’ = C: (€ — Einelastic). In the current 
approach, inelastic strains are not taken into account, i.e. Einelastic = 0. 


2.4.2 Extension with the Phase-Field Approach (Diffuse 
Cracks) 


The coupled HM formulation is extended to account for cracking with the phase- 
field modeling approach for brittle fracture. With the phase-field approach (Bourdin 
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d=0 
’ d=! 
Fig. 2.11 Discrete representation of fractures (left) and diffused phase-field representation of frac- 
tures (right) 


et al. 2000, 2008), we introduce a scalar variable that represents the damage of the 
materiald : Q — [0, 1]. A material, which is free from damage (intact) is represented 
by d = 0, while a fully degraded material by d = 1. Damage evolves in a diffuse 
manner, i.e. the crack phase-field variable evolves continuously, as shown on the 
right side of Fig. 2.11. 

We couple the evolution of the crack phase-field d with the momentum balance 
equation by multiplying a degradation function g(d) based on the evolution of d 
with the effective stresses. Before doing so, we split the effective stresses into active 
o/ and non-active o! parts computed from decomposing the elastic strain energy 
W (Amor et al. 2009; Miehe et al. 2010) as follows 


aw, = awe 
o' = g(d)—— + — 
de de 


= g(d)ol, + o! (2.7) 


Note that only the active, crack-driving part is degraded by g(d). The momentum 
balance equation is rewritten as 


0 = div (g(d)o!, +o! —apSipirl) + [1 — Ø) psr + Sidprrlg. (2.8) 

In this study, g(d) is given by 
g(d) =(1—d) +n, (2.9) 
where n is the residual stiffness (e.g., 7 = 1 x 1078 Pa) assigned in fully degraded 
part (g(d) = 0) for numerical stability. An overview regarding the derivation and 


application of a similarly combined unsaturated HM and phase-field approaches is 
described in Cajuhi et al. (2018). The crack phase-field evolution is given by 
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and n = | and 2 stand for so-called AT, and AT, models respectively (Pham et al. 
2011; Bourdin et al. 2014), and c,, = 2/3 for AT; and 1/2 for AT». 

In the phase-field literature, the split of the energy is sometimes referred to 
anisotropic response (Ambati et al. 2015). In this book, the term anisotropy refers to 
the bedding orientation and constitutive behavior. For referencing the energy decom- 
position (split), we use the term tension-compression asymmetry. There are a few 
popular split models in the literature, nevertheless, most of them can only be used for 
isotropic materials. As for anisotropic materials, the principal directions of the stress 
and strain no longer coincide. Hence, the orthogonality condition, a requirement 
for the split method to be variationally consistent, is not guaranteed. Thus, further 
treatment is required. In our recent work, we generalized two existing split models 
for isotropic materials in order to account for materials with anisotropic constitutive 
behavior. Borrowing some basic concepts from He and Shao (2019), we performed 
an energy preserving transformation so that the orthogonality condition is satisfied, 
and thus, a decomposition of an arbitrary anisotropic elastic model is feasible. We 
encourage those interested readers to follow our manuscript for a detailed clarifica- 
tion (Ziaei-Rad et al. 2023). Exemplarily, we show the verification of the framework 
in Sect. 2.5.1.2. 


2.4.3 Extension with the Discrete Element Method (Discrete 
Cracks) 


The coupled HM formulation is extended to account for cracking with the Finite 
Discrete Element Method (FDEM). 


2.4.3.1 FDEM Mechanical Solver 


Originally proposed by Munjiza (2004), the Finite Discrete Element Method (FDEM) 
utilizes principles of continuous numerical approaches in combination with those of 
discontinuous methods to simulate interaction of multiple deformable bodies (Lisjak 
et al. 2014). In other words, while the elastic deformation of discrete bodies and the 
initiation and propagation of fractures are captured using Finite Element analysis, the 
interaction between individual blocks of elements are realized by Discrete Element 
analysis. In FDEM, each block of the problem domain is discretized into a finite 
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Fig. 2.12 a FDEM discretization of the intact discrete blocks containing triangular finite elements 
and zero-thickness joint elements, b dual grid of FDEM flow solver 


element mesh with zero-thickness joint elements being inserted on its intersection 
(common) edges of the adjacent elements as shown in Fig. 2.12. Then, an explicit 
finite difference time integration scheme is applied to solve the equation of motion 
for the discretized system which reads as 


3?x ax 
M— + C— + Fin(X) — Fox (x) — F.(x) = 0, (2.12) 
dr? at 


where M and C denote the diagonal lumped mass and damping matrices, respectively. 
x is the nodal displacement vector, and Fint, Fext and F, are the nodal internal, 
external and contact forces, respectively. Contact forces, F., are calculated between 
contacting discrete blocks either in the normal direction, using a penalty method, 
or in the tangential direction, employing Coulomb-type friction law. Accordingly, 
internal resisting forces, Fint, are calculated as summation of elastic forces and joint 
elements’ bonding forces. 

Following a cohesive-zone approach, the mixed-mode crack model consists of a 
typical stress-strain curve where the bonding stresses are proportional to the sepa- 


ration and the sliding of the fracture edges. More details of the crack model can be 
found in the literature (Mahabadi et al. 2012). 


2.4.3.2 FDEM Flow Solver 


In order to account for the saturated/unsaturated porous media flow within the frame- 
work of FDEM, a vertex-centered finite volume scheme is adopted. Based on this 
approach, the conservation of the fluid phase is enforced on a dual grid of polyg- 
onal control volumes. As illustrated in Fig. 2.12, each control volume is generated 
around vertices of the mechanical solver grid connecting centers of triangular ele- 
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ments through midpoints of primary grid. The adopted mass balance equation of 
fluid phase in the porous medium is a simplified version of Eq. 2.2 with incompress- 
ible solid phase. The hydro-mechanical coupling in FDEM is realized by a staggered 
approach alternating between solid and flow solvers. Within this approach, the valua- 
tion of the stress field is affected by the pore pressure according to the effective stress 
principle, while the flow calculations are affected by the volume change associated 
with each computational node. 


2.5 Process Simulation Results 


We show different simulation setups and their results using the methods described in 
Sect. 2.4. At first, we model the mechanical processes at laboratory scale considering 
the anisotropy of the Opalinus Clay and using the FDEM formulation (Sect. 2.4.3). 
The response of the samples is compared with the laboratory experiments discussed 
in Sect. 2.3.1. Cracks are modeled either discretely or diffusely and we compare 
both approaches. At the field scale, the deformation and pore pressure behaviors are 
described with the continuum approach from Sect. 2.4.1 and the hydro-mechanical 
effects such as the shrinkage due to the de-saturation are discussed. We model the 
formation of the excavation damaged zone (EDZ) with discrete cracks (formula- 
tion from Sect. 2.4.3). Finally, desiccation cracks in the field are modeled with the 
combined poromechanical phase-field approach (formulation from Sect. 2.4.2). 


2.5.1 Laboratory Scale 


2.5.1.1 Crack Paths Due to Deformation 


As indicated by the achieved experimental results in Sect. 2.3.1.1, not only the defor- 
mation response but also the strength characteristics of Opalinus clay are strongly 
influenced by the bedded structure of the samples. In the 2D FDEM method, the 
anisotropic elastic behavior is captured by a transversely isotropic stress-strain law, 
while the strength directionality is modeled adopting a so-called smeared approach 
(Lisjak et al. 2014). Within this approach, the strength of each joint element is directly 
linked to its relative orientation with respect to the bedding layers. In this way, the 
strength parameters and fracture energies fluctuate linearly between the minimum 
limiting value, when the cohesive element is parallel to the bedding layer, and the 
maximum limiting value, when the cohesive element is perpendicular to the bedding 
layer. 

The FDEM models were calibrated to match the obtained experimental peak 
loads reported in Sect. 2.3.1.1. The laboratory-scale model along with the employed 
boundary conditions as well as the mesh details for P-, S- and Z-samples are illus- 
trated in Fig. 2.13. The experimental loading conditions were realized by explicitly 
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Fig. 2.13 a Numerical model setup of SCB experiment and mesh details for b P-sample e Z-sample 
and d S-sample 
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Table 2.3 Material and input parameters for FDEM simulations of three point bending test of OPA 


Parameter Value 
Triangular finite elements 

Density p (kg/m?) 2430 
Young’s modulus parallel to bedding E>, E3 (GPa) 10 
Young’s modulus perpendicular to bedding Eı (GPa) 4 
Poisson’s ratio parallel to bedding v23 0.35 
Poisson’s ratio perpendicular to bedding v12, 113 0.25 
Shear modulus G12, G13 (GPa) 3.6 
Crack elements 

Tensile strength parallel to bedding, ft, min (MPa) 1.3 
Tensile strength perpendicular to bedding, ft, max (MPa) 2.3 
Cohesion parallel to bedding, Cmin (MPa) 3.5 
Cohesion perpendicular to bedding, Cmax (MPa) 30 
Friction angle, ¢;(°) 22 


modeling the supports and the loading platen. In order to avoid the impact effect of 
loading platen, the loading velocity increases linearly from zero to 0.1 m/s in 1 ms. 
Following the procedure reported in the work of Tatone and Grasselli (2015), the 
calibrated mechanical material properties in the FDEM model are summarized in 


Table 2.3. 
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Fig. 2.14 Fracture path trajectory (top) and corresponding distribution of horizontal stress [Pa] 
(bottom) at 0.5 mm of crack mouth opening displacement for Opalinus clay of P-, Z- and S-samples 
with a0 = 0°, b 0 = 45° and c 0 = 90°, respectively 


Figure 2.14 illustrates the numerical simulation results concerning the fracture 
path trajectories upon failure and the distribution of horizontal stress for the three 
cases. The crack propagation shows good agreement with the experimental observa- 
tions depicted in Fig. 2.5. 


2.5.1.2 Pure Shear Test 


In order to verify the anisotropic phase-field model detailed in Ziaei-Rad et al. (2023), 
we present a set of numerical examples of a single edge notched shear test. The 
geometry and boundary conditions for this examples are shown in Fig. 2.15. We 
consider a square plate (0.001 m x 0.001 m) under plane strain conditions with a 
pre-notched crack at the middle height of the specimen. The length scale parameter 
was chosen as £ = 1 x 107 m which is small enough with respect to the specimen 
dimensions. The AT; model was used in the simulations. We used the parameters 
of Opalinus clay listed in Table 2.3. Note that the critical surface energy G. = 100 
N/m was taken to be independent of the material orientation. 

As boundary condition, we apply fixed displacement (zero in all directions) at the 
bottom edge (y = —L/2). For simulating shear, we prescribe a displacement along 
the x-direction at the top edge (y = L/2). The computation was performed with a 
monotonic displacement controlled loading with a constant displacement increment 
Au = 3 x 10°’ m. 

Figure 2.16 depicts the crack paths for orthotropic volumetric-deviatoric (Amor 
et al. 2009) split under shear loading with three material orientation angles 
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Fig. 2.15 Schematic of a cracked square plate under a displacement shear loading. An incremental 
displacement loading is applied on the top edge. The material principal axes (eı, e2) have an angle 
of œ with the global coordinate system (x, y) 


(a = —7/4,0, and x /4). Herein a significant difference in crack response is observed 
between the three cases. 

The volumetric-deviatoric decomposition assumes the volumetric compression 
does not contribute to the crack evolution. However, in case of anisotropy, the Poisson 
effect may pronounce more in a direction compared to the others, thus, even under a 
uniaxial compression, some crack driving forces can lead into damage development. 
Consequently, the crack response of a compression-dominated loading is depending 
on differences in directional stiffness of material. An implementation of the presented 
anisotropic phase-field formulation combined with the unsaturated HM approach is 
underway. 


2.5.2 Field Scale 


2.5.2.1 Hydro-Mechanical Effects and Model Response 


In this section, we numerically model the hydro-mechanical coupled effects in 
the context of the CD-A experiment using the macroscopic approach presented in 
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(a) a = -n/4 (b) a=0 (c) a= 7/4 


Fig. 2.16 Cracked square plate under shear. Phase field contours of orthotropic volumetric- 
deviatoric for three material orientations. The same parameters were used for all cases. Here, we 
observe that the crack paths are strongly dependent on the material orientation. With œ = 7 /4 there 
exist straight cracks that grow towards the lower right corner, almost along the material principal 
axis with lower stiffness. The two other cases differ from an expected isotropic response 


Sect. 2.4.1. The two-dimensional model consists of a cross section of the open and 
closed niches and surrounding domain, Fig. 2.17. We prepare an isotropic and an 
anisotropic model, in which the bedding direction is assumed to be approximately 
horizontal as monitored in the field (Ziefle et al. 2021). The OPA parameters and 
properties are extracted from several references in the literature (Bock 2000; Mar- 
tin and Lanyon 2003; Garitte et al. 2017; Bossart et al. 2018; Zhang et al. 2019) 
and presented in Table 2.4. The transversal component of an anisotropic property 
is defined with the subscript 1, while its parallel components are defined with sub- 
scripts 2 and 3. These subscripts corresponds to a local coordinate system oriented 
according to the bedding layers, i.e. 1 is orthogonal to the bedding. We compute 
the isotropic equivalent properties as the average of the anisotropic properties, e.g. 
the isotropic equivalent elastic modulus corresponds to the average of its parallel 
and perpendicular components. The isotropic equivalent properties are denoted with 
the subscript (.)iso in Table 2.4. The initial stress values are based on in-situ stress 
measurements (Martin and Lanyon 2003). Their orientation correspond to the global 
coordinate system, i.e. in two-dimensional setting their horizontal and vertical com- 
ponents op and o, lie on the global coordinate system (x, y) (Fig. 2.17), while in 
three-dimensional setting on, om and oy on (x, y, z) (Fig. 2.21). The initial pore pres- 
sure and in-situ stresses are applied on the domain. After instantaneous excavation, 
the seasonal boundary condition computed from the RH measurements (Fig. 2.1) 
with the Kelvin’s equation (Eq. 2.1) is applied at the niches. In the closed niche, 
the drops in the RH curve are related to installations and openings of the sealing 
door. We consider that the latter remains under high RH, i.e. neglect the drops, and 
assume a constant pore pressure of — 1 MPa. The open niche follows a cyclic pressure 
evolution, interpolated from the first season of the RH curve, Fig. 2.1. 

Figure 2.18 shows the saturation results using the anisotropic model for March 
2020 (first desaturation peak in winter) Fig.2.18a and September 2020 (first re- 
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Table 2.4 Parameters and material properties used for simulating the CD-A experiment 


Parameter Description Symbol Magnitude Unit 
Gravity (magnitude) g 9.806 m.s? 
Temperature To 293.15 K 
Initial in-situ stresses On —3.0 x 10° Pa 
oH —5.0 x 10° 
oy -7.0 x 10° 
In-situ pore pressure PLRO 2.0 x 10° Pa 
Grain density PSR 2520 kg-m 
Porosity $ 0.105 
Liquid density PLR 1000 kg. m~’ 
Dynamic viscosity HLR 1.0 x 103 Pa. s7! 
Bulk modulus Kı 1.0 x 10° Pa 
Intrinsic ksı 0.4 x 10-19 m? 
permeability 
ks2, ks3 1.0 x 10-19 m? 
ksiso 0.8 x 10-19 m? 
Retention curve gas entry Pe 10.0 x 10° Pa 
(van Genuchten) 
shape par. n 0.45 
res. sat. SLres 0.10 
Rel. permeability | min. Agel 1.0 x 107 m? 
(van Genuchten) 
shape par. n 0.45 
Elastic moduli Eı 6.0 x 10° Pa 
E, E3 13.8 x 10° Pa 
Eiso 11.2 x 10° Pa 
Poisson’s ratio v12, V13 0.22 
v23 0.44 
Viso 0.37 
Shear moduli G12, G13 3.2 x 10° Pa 
Biot coefficient ap 1 


saturation peak in summer) Fig.2.18b. Note that the horizontal direction is more 
permeable and, consequently, the de-saturated zone is more extended in this direction 
when compared to its orthogonal counterpart. While the saturation in the closed 
niche remains high, the saturation near the wall of the open niche increases between 
March and September. At the same time, the overall desaturated zone extends along 
the domain, but still indicates a local effect. Except for direction dependency, the 
same conclusions are drawn in the isotropic model. 

The aforementioned tendencies, i.e. de-saturation and re-saturation during winter 
and summer, respectively, have been compared with the monitored water content 
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Fig. 2.17 Setup of the two-dimensional model of the CD-A experiment consisting of open and 
closed niches. The niche radius is 1.15 m and distance between them is 10 m. The bedding direction 
is horizontal and indicated in the figure 


saturation saturation 


(a) (b) 


Fig. 2.18 Saturation changes in deformation results in the 2D model of the CD-A experiment. The 
region near the open niche dries, de-saturates, during winter and re-saturates during the summer 
months. Convergence is observed due to an increase of the deformation near the niches. Note that 
both effects are more pronounced in the open niche in comparison to the closed niche 


and electrical resistivity results. A decrease in saturation has been associated to a 
decrease of water content and an increase of the resistivity, while the opposite occurs 
when the saturation increases (Amabile et al. 2020; de Jong et al. 2020; Cajuhi et al. 
2021). These comparisons, including a geologic characterization of the niches have 
revealed that inhomogeneities like cracks, shear zones, carbonate-rich and sandy 
lenses subfacies in the claystone have a strong influence on the seasonal desaturation 
or re-saturation behavior in the near field. A detailed description of this study and its 
conclusions are reported in Ziefle et al. (2022a,b). 

The numerical results show that the seasonal saturation changes affect the defor- 
mation behavior. We compare the strain response of the isotropic and anisotropic 
models, Figs. 2.19 and 2.20, respectively, in the first winter and summer peaks. On 
both niches, the contribution of the strains due to excavation is equal. The increase 
of deformation represents convergence of the niches. Note that the projected initial 
stress is the same in both isotropic and anisotropic models. In the isotropic model, 
the strain responses are similar, while the strains in vertical direction are larger than 
in the horizontal direction in the anisotropic model. Note that the strains decrease 
from winter to summer when comparing the responses of both models. This decrease 
is specially observed in the vertical and shear components, and particularly, in the 
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(a) (b) 


(c) (d) 


Fig. 2.19 Deformation response of isotropic model. Blue regions denote negative strains, while 
red ones relate to positive strains. The strain during the first desaturation peak (March 2020) are 
shown on the left, while the strains during the second peak are shown on the right (September 2020). 
An equal strain response is observed near the niches due to excavation. The change of saturation 
induces slightly larger deformation in the vicinity of the open niche 


open niche, where seasonal effects are more pronounced and accounted for in the 
model. 

As the driving mechanism of the crack phase-field is the active elastic energy, 
which evolves with the strains (see Sect. 2.4.1), we expect to model cracks due to 
drying in the open niche using this approach. Desiccation cracks are expected to 
occur due to mode I loading, i.e. tension cracks. In the field, they were observed on 
the wall surfaces, i.e. where the rock is exposed and consequently drier in comparison 
to the undisturbed rock. Such cracks are expected to propagate from the niche walls 
into the depth of the rock. To further understand the stress distribution on the wall 
surfaces, we have prepared a 3D model of the experiment. For completeness, the 
model also considers the lock to the closed niche and the gallery that provides access 
to the niches. The seasonal boundary conditions are applied in the open niche, lock 
and gallery. The numerical model is shown in Fig. 2.21. The numerical model has 
been also utilized for interactive visualization purposes in the virtual laboratory 
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(a) (b) 


(c) (d) 


(f) 


Fig. 2.20 Deformation response of anisotropic model. Blue regions denote negative strains, while 
red ones relate to positive strains. The strains due to the first desaturation peak (March 2020) are 
shown on the left, while the strains due to the second peak are shown on the right (September 2020). 
An equal strain response is observed near the niches due to excavation. The change of saturation 
induces larger deformation in the vicinity of the open niche. Due to anisotropic elastic properties, 
the strains in vertical direction are larger than in the horizontal direction 


prototype developed by Graebling et al. (2022a) and briefly discussed in Sect. 2.6. 
As aforementioned, the bedding direction follows the horizontal direction of the 
niches’ plane. With respect to the gallery, the bedding shows an average orientation 
of (147/57) degrees (Jaeggi et al. 2020; Ziefle et al. 2021) as pointed out in Fig. 2.21. 

For a comprehensive and comparative visualization, we plot the saturation behav- 
ior in the model, Fig. 2.22. The upper part of the model is not shown in the visualiza- 
tion. Similar to the two-dimensional model, the saturation decreased in March 2020. 
A saturation decrease is observed in the open niche, lock and gallery as a result of 
the pressure changes. Due to the application of a constant pore pressure (boundary 
condition) in the closed niche, the saturation remains very high, i.e approximately 1. 
Note that the extent of the de-saturated zone in the vicinity of the gallery, lock and 
open niche remains approximately constant in size. An excavation damaged zone is 
not set explicitly in the model. 
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Fig. 2.21 Setup of the three-dimensional model of the CD-A experiment. In-situ pore pressure and 
stresses are applied as initial conditions. The closed niche is subjected to a constant pore pressure, 
while the latter varies seasonally in the open niche 
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Fig. 2.22 Saturation response in the three-dimensional model of the CD-A experiment. The satu- 
ration increases from March 2020 a to September 2020 b in the open niche and gallery 


(a) (b) 


The saturation changes and a comparison between the modeled deformation and 
monitored convergence data is found in Ziefle et al. (2022c). As also observed in the 
2D setting, the strain response of the model is affected by the excavation and by the 
seasonally change of air humidity. The strain response during the first de-saturation 
peak in March 2020 and re-saturation in September 2020 are shown in Fig. 2.23. 
Overall, all shown strain components increase up to winter and decrease in summer 
as observed, for example, in the strain component parallel to the axial direction of 
the gallery, Fig. 2.23a, b, and shear component, Fig. 2.23e, f. In the regions were the 
gallery and niches intersect, there is an increased influence of the saturation changes 
happening in the gallery. This can be visualized in all shown strain components, but 
specially in the direction parallel to the niches’ axis as depicted in Fig. 2.23c, d. 
Note that there is a development of positive strains (stretching) on the walls of the 
niches and gallery. These contribute to an increase of the energy required for crack 
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Fig. 2.23 Deformation response of the three-dimensional, anisotropic model. The results concern- 
ing the first desaturation peak (March 2020) are shown on the left, while the ones during the first 
de-saturation peak are shown on the right (September 2020). The strains are positive at the gallery 
and niches’ walls. Note that the strains are larger in the open niche 


propagation and contributes to the fact that cracks are prone to develop on the walls 
of the niches. 

The described model behavior and mechanisms serve as a basis for model exten- 
sion with respect to cracking. The proposed, extended model (Sect. ) is first 
applied in a laboratory scale test related to the desiccation of a slab, which has 
been chosen due to the similarities with CD-A. Several tests concerning the energy 
splits, model formulation and understanding the influence of the numerical param- 
eters involved in the coupled Richard mechanics phase-field approach have been 
carried out in Cajuhi et al. ( ). In the latter, the calibration and applicability of 
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the coupled framework in the field scale are described in detail. The implemented 
framework is available in the open source software OpenGeoSys (Bilke et al. 2019, 
2022). A general application of the framework in three-dimensional setting using 
a homogeneous geological window of the open niche and in context of desiccation 
cracking is discussed in Sect. 2.5.2.3. Before entering the fully coupled problem, we 
model and discuss the cracks initiated during the excavation using a FDEM approach 
in the following. 


2.5.2.2 Excavation Damaged Zone 


The Excavation Damaged Zone (EDZ) is considered as one of the main factors com- 
promising the integrity in the near field of underground repositories. In the framework 
of FDEM, simulation of EDZ formation is captured with two separate model runs 
(Lisjak et al. 2014). In the first run, the in-situ stress state is applied by translating the 
specified stresses into nodal forces. After reaching the steady state, the displacement 
of outer boundaries are fixed allowing the model to maintain the in-situ stress condi- 
tions (Mahabadi et al. 2012). The excavation step would then be initiated employing 
a core replacement technique (Kavvadas 2005). Within this approach, the Young’s 
modulus of the core material, i.e., the excavation region, is reduced in a step-wise 
manner to mimic the gradual de-confinement of the rock mass resulting from the 
tunnel advancement. 

Due to the heterogeneous geology of the Opalinus clay rock, we evaluate the 
influence of the bedding layer orientation in forming an EDZ. In the first case, 
as shown in Fig. 2.24, the numerical model setup includes a 15m x 15m square 
domain with the opening region at its center with internal radius of 1.15m. The 
geological window is representative of a one-quarter model of CD-A twins in which 
the bedding direction is horizontal. Additionally, the EDZ formation is studied with 
bedding direction inclined to tunnel axis at 45 degrees. In the second case, due to 
the asymmetric nature of the problem, a full 30m x 30m square domain has been 
considered. The employed material properties for transversely isotropic Opalinus 
clay is assumed same as the ones reported in Table 2.3. For both numerical test 
cases the vertical and horizontal in-situ stress are equal to 6.5 MPa and 4.5 MPa, 
respectively. 

Figures 2.25 and 2.26 show the gradual increase of the failure pattern at different 
ratios of deconfinement, i.e., A. Deconfinement is defined as the ratio between the 
decreased stress magnitude at the vicinity of the tunnel face to the initial magnitude 
of in-situ stress. Note that the crack propagation is parallel the bedding orientation 
in both cases. At the end of the simulation, we notice that the crack network formed 
with the horizontal orientation is more complex and connected in comparison to the 
case where the bedding is inclined with 45°. On the other hand, the fracture network 
extends deeper into the formation for the second case. The percentage of broken 
joints as a function of distance from tunnel face are plotted for both cases at à = 0.6 
in Fig. 2.27a. The steeper curve represents the 45°. Moreover, in Fig. 2.27b we show 
how the distribution of broken joints can be used in the creation of a model with 
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Fig. 2.24 Two-dimensional FDEM geometry of CD-A experiment for modeling the excavation 
damaged zone 


(a) (b) 


Fig. 2.25 FDEM results of EDZ failure pattern near the open niche with horizontal bedding layer 
at a à = 0.25, b à = 0.75 


(a) (b) 


Fig. 2.26 FDEM results of EDZ failure pattern near the open niche with 45°irc bedding layer at 
a à = 0.3, b à = 0.6 


permeability values that depend on the distance from the niche wall, i.e., the near the 
wall the more permeable the model, since more joints are broken. 

In other words, the complex crack network induced by the excavation can represent 
a zone of more instability and larger permeability in comparison to the intact rock. 
To further understand the latter, hydraulic models need to be taken into account 
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Fig. 2.27 a Distribution of broken joints percentage around underground opening. The percentage 
of broken joints in the model with horizontal bedding direction is larger (y-axis) than that in the 45° 
model. In the former, the cracks are less deep (x-axis) than those generated in the inclined model. 
Crack development motivates the EDZ abstraction concept and altered permeability profile shown 


in (b) 


and coupled with the mechanical formulation. Ongoing work is being performed to 
account for hydraulic effects in the FDEM approach. 


2.5.2.3 Desiccation Cracks 


In this subsection, we use the fundamental equations from the continuum based 
methods (Sect. 2.4) to simulate desiccation cracks in three-dimensional setting. For 
the latter, no changes in the phase-field model are required, i.e. the same evolution 
equation is utilized and there is no need to adapt the algorithm to initiate the cracks or 
due to branching. The verification of the framework is based on the setup discussed 
in Cajuhi et al. (2018). 

Since the open niche shows more pronounced hydro-mechanical effects, we pro- 
pose a quarter three-dimensional model of a so called homogeneous geological win- 
dow as shown in Fig. 2.28. The thickness of the model is the same as the radius of 
the niche (1.15 m). The total height/width is 5 m. A drying pore pressure related to 
the RH drop during the first winter is applied at the niche wall. At the remaining 
boundaries, no flow and symmetry boundary conditions are applied. Since the (free) 
movement of the niche wall, imposed by mechanical boundary conditions, can affect 
crack initiation and propagation and these conditions are not well defined in the field, 
we propose a free and a restrained setting. In the latter, the movement of the nodes 
at the front and back surfaces of the niche is hindered in the direction normal to the 
plane. 

In the following, the results concerning the degree of saturation and cracking are 
shown. For visualization purposes, a threshold filter is applied on the crack phase- 
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Fig. 2.28 Quarter model of 
the three-dimensional model 
of the CD-A experiment 
(geological window). In both 
free and restrained settings, 
symmetry boundary 
conditions are applied. In the 
restrained case, the 
movement of the niche wall 
is hindered by restraining the 
normal displacement of the 
nodes 


field variable so that the elements where damage values equal or greater than 99% 
are removed from the blended out from the mesh. When the wall is allowed to move 
freely, a single crack initiates and propagates along as shown in Fig. 2.29a. The crack 
initiates at a pore pressure corresponding to approximately 94% RH and 80-81% 
saturation. As expected, restraining the front and back surfaces, and consequently 
the movement of the wall, induces earlier crack propagation and a more complex 
crack pattern. In the latter, the crack initiates approximately at 97% RH and 90- 
92% saturation. Both cases lie in the expected drying RH range (Fig. 2.1), also 
reported in Regard et al. (2021), where crack development at the niche walls has 
been monitored and quantified. 

The current phase-field evolution equation is not able to account for the mechanical 
and hydraulic anisotropies of the material. For this reason, the developed cracks do 
not necessarily propagate parallel to the expected bedding orientation as observed 
experimentally. Ongoing work is being carried out to extend the numerical approach 
with to anisotropic setting. Basis for the extension is the framework proposed in Ziaei- 
Rad et al. (2023) and discussed in Sect. 2.4. This extension should affect the crack 
orientation and induce cracks parallel to the bedding. 

Additionally, the numerical simulation of desiccation cracking has been addressed 
in the framework of FDEM with coupled RM approach. The problem domain includes 
a 15m x 15m two-dimensional quarter model with 1.15 m niche radius. The assumed 
mechanical properties are the same as the ones considered in previous FDEM models 
dedicated to numerical simulation of fracturing in OPA. Fluid properties, porosity, 
values of permeability matrix and van Genuchten model parameters are all taken 
from Table 2.4. As mentioned in Sect. 2.4.3.2, the volume conservation equation 
and the mechanical equilibrium both are solved with finite difference explicit time 
integration scheme. This, consequently, forces smaller time steps for stability and 
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Fig. 2.29 Joint visualization of saturation and crack phase-field responses of quarter model of the 
open niche. The model desaturates near the wall of the open niche (red color corresponds to drier 
conditions, while blue color to wetter conditions). For a direct visualization, the elements with 
damage equal or greater than 99% are not plotted, i.e. the empty spaces at the wall correspond to 
cracks. A single desiccation crack initiates at the wall of the open niche in the free setting. In the 
restrained setting, a more complex crack pattern develops 
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Fig. 2.30 FDEM results of desiccation fracturing around the open niche after a 12 h and, b 24h 
simulation time 


convergence of the numerical solution. In the FDEM model, we impose a drying flux 
of q = 6.0 x 107° m/s that acts on the niche wall for 24 h. The failure pattern as 
well as the saturation distribution after 24 h is shown in Fig. 2.30. 

As can be seen, unsaturated areas, as a result of desiccation boundary condition, 
are only limited to the surrounding of the niche because of low permeability of 
Opalinus clay. Cracking starts from the top of the niche after approximately three 
hours of simulation and at approximately 3.5 MPa suction (RH between 97-98%). 
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2.6 Data Analytics and Visualization 


A prototype of the CD-A experiment using the 3D model presented in Sect. 2.5.2.1 
has been prepared and discussed in detail in Graebling et al. (2022a). The model 
is placed in the Mont Terri rock laboratory and its results can be accessed via the 
proposed interactive tool. Furthermore, experimental data and its borehole positions 
can be combined and visualized together with the numerical model. In Graebling et al. 
(2022a), convergence measurements obtained from the boreholes via extensometers 
have been compared with the strains computed with the numerical model. The direct 
visualization and comparison of the data, combined with the experience to be able 
to understand where the data come from, i.e. by visualizing the borehole positions 
and the position of the experiment itself, increases the possibilities for exchange 
between different fields of expertise and with that, a multidisciplinary approach. 
Further details on the approach and examples of the visualization are discussed in 
Chap. 5. 


2.7 Summary, Conclusions and Research Transfer 


2.7.1 Summary 


This chapter presented a workflow for understanding the Opalinus Clay charac- 
teristics, especially with respect to its shrinkage behavior. We carried out several 
laboratory experiments concerning the mechanical characterization of the rock. The 
results have shown and confirmed the strong influence of the anisotropy on the frac- 
ture behavior. In the SCB samples where the initial crack is parallel to the bedding, 
the crack propagates more easily, at a lower peak load, when compared to the exper- 
iments performed in the samples with initial crack perpendicular to the bedding. 
The tests have also shown the combined interaction of bedding and heterogeneities 
in the OPA since the cracks do not propagate uniformly in all samples. We also 
computed these experiments using FDEM and phase-field models coupled with the 
mechanics. Experience has also been gathered with respect to sample preparation 
and cutting procedures. For increasing the understanding on the effect of humid- 
ity variations, tests under different humidity conditions have been carried out. The 
tests were performed in a customized desiccator cell which has been equipped with 
different sensors to measure the controlled RH and temperature. In summary, the 
experimental tests show that an increase of suction (decrease of RH) leads to an 
increase of the strength of the material—as long as no fractures are induced that lead 
to an overall weakening of the material. The obtained data can be used as input in the 
simulations of drying/wetting and will affect the fracture behavior of the material in 
coupled hydro-mechanical simulations. 

The laboratory tests have a direct connection with the field, since the samples 
were obtained from drilled boreholes in the vicinity of the in-situ CD-A experiment. 
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In this region, the main lithology is characterized by the sandy facies of the OPA. 
We presented and used material properties obtained from the literature and from 
the experiments performed in laboratory. These data have been used as input in the 
various presented numerical models. The monitored data at the CD-A twin niches 
have been used for comparison with water content and convergence data. These 
results have been described briefly and are object of further investigations (Ziefle 
et al. 2022b, c). The twin niches at the CD-A experiment provide a good platform for 
comparing the hydro-mechanical effects under seasonal (open niche) and controlled 
high RH conditions (closed niche). Furthermore, the niches are equipped with sev- 
eral monitoring sensors, e.g. extensometers, permeability, water content, crackmeter, 
among others. The open niche shows more pronounced shrinkage effects in compar- 
ison to the closed niche, for example with respect to desiccation crack development. 

The development of excavated damaged zones and further hydro-mechanical 
effects near the excavations, for example the initiation and propagation of desic- 
cation cracks during ventilation, are motivating facts for developing and extending 
numerical approaches with respect to cracking and evaluating their applicability in 
the in-situ scale. With this objective in mind, we proposed continuum and discontin- 
uum based models. A comparative study between these approaches for the laboratory 
scale is shown in Sattari et al. (2022). We described a continuum approach for mod- 
eling the unsaturated HM-coupled problem that has been implemented in Bilke et al. 
(2019, 2022). This approach has been applied in two- and three-dimensional set- 
tings. Aim of the simulations was to understand the effects of de-saturation and 
re-saturation in isotropic and anisotropic conditions for the hydraulic (permeability) 
and the mechanic (elasticity properties). The isotropic properties have been com- 
puted as an average of the anisotropic ones. We have shown the results for the 
first de-saturation peak, i.e. largest suction, and for the first re-saturation peak, i.e. 
largest pore pressure. The de-saturated zone increases up to the smallest pressure 
and part of it remains up to the re-saturation peak. We have observed the influence 
of the saturation changes on the deformation behavior and have stated that these are 
more pronounced in the open niche, where the seasonal behavior plays an important 
role and desaturation is stronger. These conclusions are valid in both isotropic and 
anisotropic models. However, in the latter, the desaturated zone extends further due 
to larger permeability along the horizontal direction. 

The numerical results concerning the hydro-mechanical effects in 2D setting have 
shown that the contour and vicinity of the open niche are drier in comparison to the 
domain as a whole. This represents a driving mechanism for desiccation cracks. To 
obtain further information on the drying behavior of the walls, the model has been 
extended to 3D. This model is composed of open niche, closed niche, lock to the 
closed niche and gallery. We have observed an increase of the strains on the walls 
of the niches and gallery, with more pronounced effects in the latter and in the open 
niche. 

Note that the models discussed up to now have not considered the development and 
extent of an EDZ. For modeling the crack networks formed during the excavation, we 
have used an FDEM approach. The approach combines the finite element analysis 
with a discontinuous representation of the crack, i.e. through zero-thickness joint 
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elements that are inserted in the mesh. We observed a more complex crack network 
when modeling an horizontal bedding layer on the cross section of the niches that 
might represent a zone of weaker stability and larger permeability. 

The unsaturated HM approach has been extended to account for cracks due to dry- 
ing. The model has been applied into an homogeneous geological window, consisting 
of a quarter of the open niche. The free and restrained niche wall movement have 
been simulated to account for the full range of boundary conditions in the field. In 
the former, a single crack has been formed, while the latter induced a more complex 
crack pattern. The crack propagation is not direction dependency The effect of the 
hydraulic and mechanical anisotropies have not been taken into account in this model. 
The obtained crack patterns do not follow the bedding direction since the model is 
isotropic. Nevertheless, the quarter model has been able to predict desiccation cracks 
in the expected and monitored RH range (94-97%). 

Finally, an interactive visualization of data in the field proportions has been pro- 
posed. The 3D model presented in this chapter has been inserted in this virtual 
prototype of the Mont Terri Laboratory, where the CD-A experiment takes place. 
The numerical results of saturation and deformation can be rendered and directly 
compared with borehole data, its position and monitored data provided by the sen- 
sors. 


2.7.2 Conclusions and Research Transfer 


We have proposed a workflow consisting of laboratory, field scale experiment and 
visualization of data. The laboratory experiments have provided new data related to 
the mechanical and fracture behaviors of the sandy facies of the OPA. Furthermore, 
a relationship between peak load and suction has been established for this material. 
The data obtained experimentally can be used as input for coupled hydro-mechanical 
simulations. For that, e.g. the strength needs to be determined for each tested suction 
range. An important model improvement would be to use strength-suction dependent 
material properties for simulating the material in laboratory and field scales. 

The FDEM method showed good agreement with the measured laboratory experi- 
ments and has been also tested in the field scale for dealing with the EDZ. The numer- 
ically obtained EDZ shows complex crack networks, which might have an important 
effect on the permeability changes near the excavations. For further understanding 
this behavior, the hydraulic effects need to be coupled with the fracture-mechanical 
approach. This is object of ongoing work. A comparison between the response of the 
coupled hydro-mechanical methods (continuum and discontinuum based) is under- 
way (Sattari et al. 2022). With that, we will evaluate the computational costs related 
to the inclusion of the discrete elements (joints) and of the addition of the phase-field 
evolution equation in each formulation. 

The extension of the continuum framework based on the unsaturated HM approach 
using the phase-field model for brittle fracturing has the advantage and ability of sim- 
ulating desiccation cracks. For a more resolved crack, the meshes need to be refined 
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and this fact is directly connected to increased computational effort. Methods for 
automatic mesh refinement need to be further developed and explored in the context 
of HM coupled simulations. In the context of the unsaturated HM approach, the use 
of the presented anisotropic phase-field formulation is being evaluated and imple- 
mented in the open source framework. Furthermore, ongoing research is planned to 
account for swelling effects and crack closing as well as cracks due to swelling. 

The long-term planning and execution of the CD-A experiment is an important 
aspect for testing the developed/extended modeling approaches and comparing its 
results to field data. Furthermore, the interactive data visualization of in-situ data is 
an important tool for communication and data exploration. The prototype presented 
in Graebling et al. (2022a) can be extended to account for further boreholes and 
sensors. 

An interesting aspect for further research is the application of the presented 
methodology on other types of clays and rocks. This methodology can aid on finding 
the links between the different spatial and temporal scales. The combination of labo- 
ratory and field experiments, fundamental and applied research as well as interactive 
visualization represents an added value for different applications such as those in the 
geo-energy and radioactive waste disposal. 
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Chapter 3 A) 
Pathways Through Pressure-Driven de 
Percolation in Salt Rock 


Markus Knauth 


3.1 Hydro-Mechanical Properties of Salt Rocks and the 
Concept of Pressure-Driven Percolation 


Under undisturbed conditions and in the absence of significant impurities or inter- 
calations (e.g. anhydrite) rock salts are considered to be tight against fluid or gas 
pressures, which stems from their crystalline structure composed of strongly vis- 
coplastic, impermeable salt grains grown together at their grain boundaries (Mink- 
ley et al. 2013, 2015). Thus, its small natural porosity of 0.1-1% consists mainly 
of isolated brine pockets as evidenced e.g. by microscopic investigations of thin salt 
slivers. There are however stress- and pressure-conditions under which this hydrauli- 
cally inactive flow network of connected grain boundaries can be opened. Then—and 
only then—a significant increase in macroscopic permeability can be observed due 
the creation of a connected flow path by opening the grain boundaries for permeation 
(Fig. 3.1). The conditions for such an increase in permeability can be classified into 
two principle mechanisms: 


e damage induced creation of internal fractures by deviatoric loading 
e pressure-driven opening of grain-boundaries by an applied fluid/gas pressure 


The first item describes the creation of flow paths by damaging the salt due to loading 
above its threshold for damage initiation (the so-called dilatancy boundary). The sec- 
ond mechanism—which is the main focal point of this study—describes the ability 
of fluids and gases to “open” the initially impermeable grain boundaries when the 
applied pressure becomes higher than the normal stresses acting on them. The attack- 
ing pressure thereby allows the medium to push its way into the grain boundaries 
and creates a connected flow path structure. Since this formation of a flow network 
occurs only after exceeding a certain threshold related to stress and pressure, this 
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Fig. 3.1 Percolation of a colored solution on grain boundaries in hydro-mechanically coupled 
triaxial compression tests on salt specimen 


effect has been termed “pressure-driven-percolation” in the field of salt mechanics. 
In the assessment of barrier integrity in salt mining—i.e. the preservation of afore- 
mentioned sealing capacity of the salt units—these two mechanisms have lead to the 
development of two associated criteria for their evaluation (Kock et al. 2012): 


e Dilatancy criterion 
e Minimum principle stress criterion (MPS) 


The dilatancy criterion describes the amount of mechanical damage the salt material 
has experienced, while the MSP compares the acting minimum principle stress with 
a theoretically attacking fluid pressure, e.g. a column of groundwater above the 
salt top or the gas pressure within a salt cavern. While the dilatancy criterion is 
typically more relevant for the assessment of damage at the contour of underground 
drifts and caverns (the excavation disturbed zone—EDZ), the MSP generally has 
a more far-reaching impact since the large-scale stress redistributions around a salt 
mine may reach far into the salt, potentially threatening its capacity as a hydraulic 
protection layer. It is important to note, that the minimum principle stress criterion is 
not the result of ahydro-mechanical modeling, but solely a stress- and pressure-based 
assessment of a given stress-state from a geomechanical modeling. As such it must 
be ensured, that it is applied at the most unfavorable state of the modeled system. 

Additionally, the MSP is conservative in the sense that it only considers the 
scalar quantity of the least compressive stress. In reality, numerous laboratory hydro- 
mechanically coupled percolation studies have shown that the pressure-driven per- 
colation is a highly oriented process, where the fluid generally percolates in a plane 
perpendicular to the least compressive stress (Kamlot 2009) (Fig. 3.2). Therefore 
it is possible that the MSP becomes too conservative, when the minimum principle 
stress may be lower than the attacking fluid pressure, but the direction of the stress 
field is such that this does not lead to a loss of the hydraulic protection capacity of 
the salt layer. 
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Fig. 3.2 Highly anisotropic percolation in salt rocks perpendicular to the least compressive stress 
under true triaxial compression (Kamlot 2009) 


Based on the previous considerations, the idea suggests itself to model the 
pressure-driven percolation as a diffusion process using a dynamic, full permeability 
tensor based on the fluid pressures and the stress tensor. However, this seemingly 
simple idea of modeling highly anisotropic diffusion leads to a number of numerical 
challenges in this and also in other scientific fields. The subsequent section will there- 
fore outline the principle idea and the arising problems more detailed, before propos- 
ing a mesh-free numerical approach taken from the field of magneto-hydrodynamics 
in astrophysics. 


3.2 Development of a Mesh-Free Percolation Model 
Inspired by Magneto-Hydrodynamics 


3.2.1 Numerical Challenges in Modeling Highly Anisotropic 
Permeation 


Most geomechanical codes capable of modeling a Darcy-type porous diffusion use 
finite-volume techniques (Moukalled et al. 2016), where the fluid resides in control 
volumes and each of these volumes may exchange fluid with its neighbors. By solving 
a balance equation and, the flow rates between those elements are determined and 
the corresponding changes in saturation and/or pressure are calculated. In order to 
determine the pressure gradient between to connected volumes i and j, a so-called 
two-point flux approximation (TPFA) is made (Aavatsmark 2002): 


Pj Pi 


Vpij = as 
ij 


€j 
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This widely-used approximation is one of the main reason for excessive numer- 
ical (i.e. false) diffusion when a highly anisotropic permeability tensor is used. 
For anisotropic permeability tensors, this approximation only yields accurate results 
when the directions of the element connection vectors (e;;) are aligned with the Eigen- 
vectors of that particular permeability tensor, i.e. its principle directions. Assuming 
that the permeability tensor is not constant and may change its principle directions (as 
it is proposed here by coupling it to the fluid pressure and stress tensor) this condition 
will never be met in the simulations. The introduced numerical error is of order n(0), 
meaning it is inherently tied to the model structure and does not converge away with 
increasing mesh resolution. It leads to numerical instabilities and implausible results 
including: 


e violation of non-negativity (may yield negative temperatures in numerical model- 
ings of anisotropic temperature distributions (Sharma and Hammett 2007)) 

e violation of the maximum principle (in a closed system, the method may yield e.g. 
higher fluid pressures than the ones it started with) (Gao and Wu 2013) 


Even though this numerical problem is known for decades and the literature shows 
numerous correction approaches in a wide spectrum of applications from thermal 
propagation to astrophysics, it is still the default approximation in many—if not all— 
established geomechanical codes. In the context of geomechanical barrier integrity 
not correcting for this problem leads to uncontrolled numerical diffusion, since the 
anisotropic fluid distribution will not stay stable and inevitably fill the whole model 
domain, thus rendering the approach useless for the evaluation of a hydraulic barrier. 


3.2.2 Basic Approach and Formulation 


In order to discretize the problem, a choice has to be made on how to discretize the 
model domain volume into a set of cells or particles i with coordinates x;. Instead 
of using a typical mesh e.g. of interconnected hexahedra or Voronoi elements, the 
approach formulated by Hopkins et al. (Hopkins 2015, 2016) considers a mesh-free 
alternative. A differential volume d*x at coordinate x. This differential volume can 
then be partitioned fractionally among its nearest particles by using a weighting 
function W depending on the distance x — x; and a “kernel size” h(x), thereby 
associating a fraction W;(x) of the volume d*x with particle i: 


1 
Wi(x) = FC — xi, h(x)) 


olx) = J (x - xi, h(x) 
j 


W(x — x;, h(x)) can be any continuous, symmetric and compactly supported func- 
tion. Its absolute value is irrelevant due to the normalization. This type of model 
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discretization therefore conceptually sits somewhere between a Voronoi tessellation 
and smooth-particle hydrodynamics. Hopkins et al. (Hopkins 2015) show, that this 
partition of unity discretization makes it possible to formulate the balance equation 
of property U for a given particle i with associated volume V; as a Godunov-type 
finite-volume equation: 


d 
VUN) + Fi Ay = 0 
j 


with the “fluxes” Fj; in or out of an effective face area A;j. 

Another crucial advantage of this method is the fact that that type of discretiza- 
tion easily allows for a second-order accurate approximation of the gradient of any 
property f as a least-squares fit over all points in the kernel region, which effectively 
makes it a type of multi-point flux approximation method (MFPA): 


VAT = F -AWE æj- xD? oid 
j 
Wir =} aj- xE- xoja) 
j 


Is it important to note, that this gradient estimator is second-order accurate for an 
arbitrary particle configuration and therefore particularly suited for the mesh-free 
method. While the method proposed in Hopkins (2015) is constructed for very general 
hydrodynamics applications and leads to solving the associated Riemann problem, 
we now differ from the original formulation by transferring the general notion to 
the field of diffusion in porous media. In a simple explicit time-stepping scheme, 
we calculate the fluxes F;; between each particle and its associated “neighbors” by 
evaluating the approximated gradient at the midway point between each particle and 
summing these contributions. 


1 
Fij = Ti -Vpij) > Aij 


In order to stabilize the anisotropic diffusion equation, this flux is compared to the 


direct two-point flux 
aoe 1 Pj — Di 
Fa ct K J l A 
ij H | i ( dij ey) i] 


Using this direct flux, the final flux F, j is given by 


Bo 0 if sgn( Fire 5 F;j) and | 29% | Fi;| 
2 F;; otherwise 


52 M. Knauth 


Additionally, a basic saturation concept is introduced, giving each particle a porosity 
value and ensuring, that no flow can occur out of an unsaturated zone by multiplying 
the calculated flow rate with an empirical saturation function as s; j used in the 
geomechanical modeling codes of Itasca (Itasca 2022a, b): 


Sij = s? . (3 — 2s) 


where s is the saturation of the zone where the fluid is coming from (upwind scheme). 
Using this approach the volumetric flow rate Q; for particle i is 


Qi=) Fj si 
J 


If the saturation at particle i is less than 1, this flow rate is used to update the saturation 
by 
Qi 


Nie Vi 


s(t + At) = si(t) + - At 


with porosity n; and associated volume V;. If the particle is already fully saturated, 
the associated pressure p; is updated instead using the fluid bulk modulus B finia: 


fluid * Qi 


B 
pit + At) = p(t) + - At 
nj - Vi 


This algorithm was implemented as a Python module in FLAC3D, since this is the 
framework exposing the relevant internal variables of the corresponding numerical 
FLAC3D model in order to have access to the corresponding stresses and strains 
from the mechanical calculations of the salt behavior. In a staggered approach, the 
time-dependent mechanical problem is solved iteratively by FLAC3D for a given 
amount of time before the new stress tensor field is then given to the newly devel- 
oped Python module for the anisotropic fluid flow calculation for the same timespan 
(albeit obviously using a different timestep for the fluid flow calculations). Although 
it is has experienced some optimizations, the Python implementation is—by nature 
of the interpreted programming language—much slower than a comparable C++ 
implementation using HPC-ready sparse matrix library would be. This would how- 
ever need to be encapsulated again as a Python module in order to continue to work 
in conjunction with the FLAC3D software. Such efforts were tested in principle with 
the framework of this project and are (in addition to optimizations of the numerical 
structure itself) a potential field of improvement for future work. Together with the 
dynamic calculation of the local permeability tensor based on the stress tensor and 
the fluid pressure (see Sect. 3.3.3), we will subsequently refer to the new Mesh-Free 
Percolation technique as the MFP-approach. 
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3.3 Validation in Analytical Test Cases and Idealized 
Percolation Situations 


In order to demonstrate the capabilities ofthe new method and also show the unphys- 
ical numerical problems arising in standard implementations, a number of test cases 
have been assembled as a benchmark for the modeling of highly anisotropic porous 
flow. These examples and their comparative results are discussed hereafter and may 
aid in validation and verification for similar methods in the context of geomechanical 
percolation processes. 


3.3.1 Diffusion of a ID-Step Function in Different Constant 
Permeability Fields 


In order to validate the basic diffusion capabilities, a 1D-step function is initialized 
in a constant isotropic permeability field (3.3) with a pressure p = 1 MPa on the left 
half and zero pressure on the right half of the model. The mesh has been created delib- 
erately in a tilted fashion (with respect to the orientation of the quadratic elements) 
in order to create the most unfavorable condition for the subsequent analysis with an 
anisotropic permeability. For the isotropic permeability k in the first modeling, the 
pressure step will diffuse and soften the transition with time given by the analytical 
solution for the pressure profile p,,, for a given time t: 


p p X — Xo 
‚r>0 Erf 
PGE eS tg [= | 


with 
— Bria: K 


nu 


where n is the rock porosity and u the viscosity of the flowing medium. Figure 3.3 
shows the pressure profile for a time of t = 600 s using the new MFP-approach in 
comparison with the native FLAC3D fluid logic. As to be expected, both codes are 
in good agreement with the analytical solution for the isotropic case. 

In the second step, we now assign a anisotropic permeability, which is only 
nonzero in the vertical direction (kų = Om?). Therefore, the numerical simulation 
should keep the initial step function perfectly intact, since there should be no flow 
into the x-direction. Again we examine the pressure profiles for both codes in Fig. 3.4 
for a time oft = 12000 s. It can be clearly seen, that the MFP-approach coincides 
with the analytical (=initial) pressure profile, while the FLAC3D model exhibits the 
numerical errors discussed in the previous sections: The fluid front does not stay 
stable and move into the x-direction while also showing higher pressures than the 
initial condition on the left side and negative pressures on the right hand side. This 
is a clear example of the undesired effects using lower-order flow approximations 
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1D Isotropic Diffusing Step Function (t=600s) 
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Fig. 3.3 Model setup and mesh for the 1D isotropic diffusing step validation example (top) and the 


modeling results of the resulting pressure profile for different approaches (bottom) showing good 
agreement with the analytical solution 


in highly anisotropic settings and is therefore proposed as a first and simple validity 
check for anisotropic approaches on porous flow. 


3.3.2 The Diffusing Ring Problem 


Another numerically challenging example for highly anisotropic flow is the so-called 
diffusing ring. Following the description in Hopkins et al., we initialize a purely 
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1D Anisotropic (Non-)Diffusion of a Step Function (t=12000s) 
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Fig. 3.4 Resulting pressure profiles for the 1D anisotropic diffusing step function, which should 
stay perfectly immobile. The MHD approach successfully captures this, while the standard two- 
point flux-approximation implemented in FLAC3D fails to keep the fluid front stable and also 
exhibits numerical problems 


azimuthal pressure field in a 2D cubic box. With cylindrical coordinates centered on 
the origin we initialize p(t = 0) = po(exp—(1/2)[(r — ro)-/ör} + &°/863]) with 
p0 =,ro = and do =. Since the permeability field is purely azimuthal, the pressure 
should then only diffuse along a ring-shape around the center. For short times, this 
problem has an exact solution: 


p(t > 0) = po(exp —(1/2)L(r — ro)” /ôr + 67/5¢7]) 


with 
86° (r, t) = 565 + 2Kr "rt 


At later times, the diffusion from both sides of the ring intersects and there is no 
longer an exact solution. Again we test this problem setup using the MFP-approach 
and the native FLAC3D implementation. As an example of the modeling results, 
Fig. 3.5 shows the pressure distribution along the annulus after a short time (t = 
60 s) in comparison to the original distribution. Even for this very short time, it 
can be seen that the MFP-results are again in good agreement with the analytical 
solution, while the FLAC3D results is already starting to deviate noticeably. For 
longer times, the pressure distribution wraps around as to be expected and stays fairly 
stable in the MFP modeling, while it again exhibits strong numerical diffusion in the 
FLAC3D implementation, eventually filling and equalizing the pressure in the whole 
domain. Therefore, this second—more elaborate—numerical example has confirmed 
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Fig. 3.5 Resulting pressure profiles for the 1D anisotropic diffusing step function, which should 
stay perfectly immobile. The MHD approach successfully captures this, while the standard two- 
point flux-approximation implemented in FLAC3D fails to keep the fluid front stable and also 
exhibits numerical problems 


the shortcomings of the default methods and the capabilities of the MFP-approach 
to accurately model this challenging setup. 


3.3.3 Inclusion of Pressure- and Stress-Dependent 
Anisotropic Diffusion 


So far, the presented results were carried out for constant permeability problems 
with different boundary conditions and setups with the intent to verify the basic 
capabilities of the proposed method to accurately model anisotropic fluid flow. At 
this point we now include the mechanisms described in the introductory section in 
order to model the pressure-driven percolation in salt rocks. 

In its principle axes, we formulate the anisotropic diffusivity tensor K in the 
following way: 


ka + k3 0 (0) 
K= 0 k+k 0 
0 0 kı + ko 


k= A(p+o;)forp>o; 
"7 [0 otherwise 


where o; are the principle stresses and p is the fluid pressure. By this construction, 
e.g. a pressure p exceeding the minimum principle stress o3 will lead to an increase 
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< FLAC3D native x MHD Implementation 


Fig. 3.6 45° degrees inclined pressure-driven percolation modeling with strong numerical diffusion 
and invalid pressures (holes) in the default FLAC3D implementation and stable percolation in the 
novel MHD approach 


in permeability in the directions of o; and o2, i.e. the fluid may move in a plane 
perpendicular to the principle stress 03. 

Using this straightforward approach we can then test its functionality by modeling 
fluid injection tests in a cubic sample subjected to different stress conditions similar 
to the ones described in Kamlot (2009). On the purely cartesian mesh (20 x 20 x 20 
elements), we initialize a fluid pressure of p in the center of the model, whose stress 
state is given by o1 = o2 = 11 MPa and o3 = 9 MPa. The given stress field was again 
deliberately chosen to lead to an unfavorably inclined fluid flow forming an angle of 
45° against the principle axes of the mesh. Since the initial fluid pressure is below 
the minimum principle stress in the sample, there should be no flow occurring until 
the fluid pressure—which is increased stepwise—reaches and exceeds that mininum 
principle stress. 

As expected, the initialized fluid stays stable until it exceeds the minimum prin- 
ciple stress. Upon exceeding the percolation front develops in the expected inclined 
orientation and is—most importantly—stable with respect to potential erroneous 
numerical diffusion in the direction of 03 (Fig. 3.6). Therefore, the proposed method 
is capable of both modeling both constant anisotropies as well as with the additional 
interaction of a spatially varying stress- and pressure-dependent permeability tensor. 


3.3.4 Application on a Large-Scale in Situ Borehole 
Pressurization Test in a Salt Mine 


A numerical recalculation of a large-scale in situ test was to be performed in order 
to demonstrate the capabilities and accuracy of the new approach. In this in-situ 
test, a 60 m long borehole with a diameter of 1,35 m was drilled between to mining 
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Fig. 3.7 Location of the test site in the mining field Springen (left) and construction of the large- 
scale 60 m borehole with 1.35 m diameter between mining levels (right) 


levels of a salt mine and additionally monitored by a seismic array (Fig. 3.7). The 
borehole was cemented in the lower part, leaving a 40 m high cylindrical cavity 
which was then pressurized with gas. In accordance with the expected stress field, 
the salt borehole was gas-tight until the gas pressures reached a level above the 
minimum principle stress at the lower end of the borehole. The thereby created 
percolation front started in the bottom part of the borehole (at the cement plug) and 
progressed in a nearly horizontal plane above the underground mine drifts (Fig. 3.8), 
which was to be replicated by the newly developed method. The hydro-mechanically 
coupled calculations of the percolation reaction in response to the stepwise increase 
of borehole pressure successfully confirm the in situ observations. Using the novel 
MHD-approach, the modeled borehole is tight against the attacking gas pressure until 
the minimum principle stress is exceeded, which happens first in the lowest part ofthe 
borehole (Fig. 3.9). The direction of the percolation plane is oriented perpendicular 
to the smallest compressive stress and therefore initially nearly horizontally aligned. 
The percolation front in the in situ—test eventually hit a almost perfectly horizontal 
weakness plane and is therefore slightly more localized, which is not included in 
this modeling. However, the new method has shown to be well suited to model the 
pressure- and stress-tensor-dependent salt percolation process both in laboratory and 
in situ scale remarkably well. 


3 Pathways Through Pressure-Driven Percolation in Salt Rock 59 


Cluster 


Pressure (bar) 
& 
° 


> 025. 
1* floor: 275 m O wrge diameter 
borehole (to-scale) 
220; D montonng borehete 
i 23.01.2012 amoreod) 
© location of seisme 
overt 
os migration front 
4 extension 
810) 
Lo open E 
> shaft Z 
are — >05 
s0 
sealed 800 
shaft 
w o 
2% floor: 360 m | 
™ A a 
kJ o ro w » 100 ix 190 
klini a 6 7 ” e 85 s ss 


Xin m 


Fig. 3.8 Pressure steps of the in situ—test (top left) and AE-measurement (setup top right) of the 
percolation front after the 68 bar pressure step in side (bottom left) and top view (bottom right) 
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Fig. 3.9 Modeled percolation front starting to develop only after exceeding the necessary perco- 
lation pressure and permeating perpendicular to the least compressive stress 
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Chapter 4 ®) 
Pathways and Interfaces Under Stress de 
Redistribution 


Markus Barsch, Thomas Nagel, Holger Steeb, Patrick Schmidt, 
Dongwon Lee, Carlos Guevara Morel, and Jobst Maßmann 


There are several finite element-based modeling approaches to deal with fissures, 
fractures and discontinuities in rocks. A recent overview can found e.g in Mollaali 
et al. (2022). As one of them, lower-dimensional interface elements (LIE) were 
implemented in OpenGeoSys! previously (Watanabe et al. 2012; Yoshioka et al. 
2019) to map fractures and fissures discretely. 

For example, in a 3D model consisting of a matrix and a fracture, the matrix 
is represented by 3D elements and the fracture by 2D elements. Crack growth can 
occur along pre-defined interfaces, i.e. the 2D elements, and therefore, the direction 
of crack path is predefined in contrast to other approaches resting on phase fields, 
extended finite elements or remeshing strategies. 


‘https://www.opengeosys.org. 
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Further, lower-dimensional approaches taking into account thehydromechanics of 
deformable fluid-saturated fractures in porous media have been implemented in open- 
source software packages like DUNE? or FEniCS? (Schmidt and Steeb 2019; Schmidt 
et al. 2021). Numerically efficient coupling approaches for fractures and the sur- 
rounding porous rock allowing for HPC computations based on lower-dimensional 
models have been successfully applied to complex 2D/3D problems (Schmidt et al. 
2022b) using modern code-coupling tools like preCICE.* 

This general framework can be applied to a range of problems aside from frac- 
ture mechanics, including fracture slip problems and geotechnical problems with 
frictional interfaces, such as soil-structure interaction. In the sequel, we briefly sum- 
marize the theory and give a short insight into selected application examples from the 
GeomInt2 project. In cooperation with researchers from STIMTEC (Renner 2020) 
and others, those simulation tools have been used for the numerical interpretation 
of pressure transients obtained from large-scale experiments in the field, like e.g. 
harmonic pumping tests at the underground laboratory Reiche Zeche, Freiberg, Ger- 
many (Schmidt et al. 2021) or at the Grimsel Test Site, Switzerland, (Schmidt et al. 
2022a). Based on those numerical-experimental investigations it was shown that 
fully-coupled hydromechanical approaches give a more physically sound insight into 
pressure transients especially if results are compared to simpler and more “classical” 
diffusion-based modeling approaches. 


4.1 Hydro-Mechanics of Deformable High-Aspect Ratio 
Fractures 


Viscous fluid flow in fractured porous rocks is causing strongly coupled hydro- 
mechanical interaction phenomena by local aperture, i.e. volume, changes of the 
fractures and resulting variations of the fluid pressure state. These phenomena have 
been investigated by various researchers over the last decades. For an overview we 
refer to a recent PhD thesis written during the period of the GeomInt and GeomInt2 
project (Schmidt 2022). 

Modelling approaches of fractured reservoirs require an efficient description of the 
hydro-mechanical processes. This includes on the one hand the accurate description 
of the deformation of the fracture and on the other hand the dynamic exchange of 
momentum between the fractures and the surrounding rock matrix. From a continuum 
perspective, model assumptions have to be made for the geometry of the fracture, 
and in a lower-dimensional setup, for the velocity profile within the fracture domain. 
For low Reynolds-number flow (Re < 1), a hybrid-dimensional model approach is 
obtained by an integration of the assumed parabolic flow velocity profile in normal 
fracture direction €3. This “a priori” integration is reducing the fracture flow domain 


? https://www.dune-project.org. 
3 https://fenicsproject.org/. 
4 https://precice.org. 
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Br 


Fig. 4.1 Representation of a single fracture Tr, embedded in a poro-elastic material body B?¢ 
with material surface r?°. For the kinematical description of the fracture, a local basis system 
€; is introduced. The basis vector €3 is aligned with the fracture surface normal. The mechanical 
deformation of the fracture, and therefore the local volume change, is depicted with the evolution 
of the aperture ô, cf. Schmidt et al. (2022b), Schmidt (2022) 


by one dimension, cf. Fig. 4.1. The resulting model is derived from 3D continuum 
mechanics by a consistent evaluation of the conservation of mass and balance of 
momentum under fully saturated conditions (Vinci et al. 2014; Schmidt 2022). 
Without going into details of the derivations of models, we recapitulate here only 
the main governing equations. The bulk properties of the porous rock are described 
with the quasi-static linear poro-elastic equations (Biot 1941; Steeb and Renner 2019) 


—divo = pb, 
lap kf du, 
— P _ Z div gradp = —a div —, a) 
M at vi ðt 


In Eq. 4.1 we introduced the local storativity 1/M of the porous rock, Biot’s coupling 
parameter œ, the isotropic hydraulic conductivity k/, and the effective weight of the 
fluid y£ ® The total stresses are given by ø, the pore pressure is denoted with p, body 
forces are introduced as p b and the displacements of the solid skeleton are given by 
u,. The set of coupled equations are supplemented with boundary conditions (Steeb 
and Renner 2019). 

From first principles, a pressure diffusion equation could be derived for the evo- 
lution of the pressure state within the fracture. For high aspect ratio fractures, cf. 
Vinci et al. (2014), the final partial differential equation could be written in a simple 
non-linear format taking into account the evolution of aperture ô as well as fluid leak- 
off through an exchange term with the matrix g (normal component of the seepage 
velocity) 


ap 8? 1098 å 
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In Eq. 4.2, we have introduced fluid properties like the effective dynamic viscosity 
n/R, the inverse of the fluid’s bulk modulus £^ and the fluid pressure in the fracture 
p. Details concerning the derivation and the mentioned dimensional analysis of Eq. 
4.2 could be obtained from the references mentioned e.g. in Schmidt (2022). The set 
of coupled partial differential equations in weak form for the poro-elastic domain 
and the embedded fractures could be implemented into a Finite Element framework, 
e.g. within software packages mentioned in the introduction. 


4.1.1 Stiffness of Fractures 


The change of deformation of a fracture does not only depend on the fluid pressure 
transients p(X, t). The local deformation of the fracture is also affected by local 
contact forces between asperities of the solid surfaces. In case of natural fractures 
which have always rough fracture surfaces, those contact forces could for instance 
occur in the case of “mechanically closed” but still “hydraulically open” fractures. 
In the present lower-dimensional modelling approach, the corresponding effective 
fracture stiffness caused by local contact forces could be taken into account and 
expressed with the contact related normal stresses (Bandis et al. 1983) 


op = EF ———— 0. (4.3) 


Note, that in Eq. 4.3 we considered only normal stresses, i.e. stress components 
normal to the fracture plane. Shear stresses within the fracture surface are neglected 
a priori. In the most simple linear constitutive approach, the effective contact normal 
stresses are proportional to the change of aperture described with the fracture closure 
6° and the maximum fracture closure ôhax» cf. Gens et al. (1990) and results in Fig. 
4.2 . In this linear constitutive relation we have introduced, as an inherent material 
parameter of the contact law, an effective normal stiffness E™. 

On the laboratory scale, the introduced effective stiffness parameter could be 
determined in well-defined experiments. Therefore, natural fractures, obtained from 
mode 1 hydraulic fracture experiments, are mechanically tested, e.g. under harmonic 
excitation, cf. Fig. 4.3. 


4.1.2 Strong Versus Weak Coupling 


Implementation of the coupled set of partial differential equations introduced in 
terms of the biphasic poro-elastic and hybrid-dimensional flow model (Eqs. 4.1 
and 4.2) might be distinguished regarding the chosen coupling scheme. In terms 
of the numerical strategy to reach equilibrium, chosen numerical coupling schemes 
are distinguished, cf. Fig. 4.4. Here, we may distinguish between staggered and 
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Fig. 4.2 Left: Relative normal stress changes as a function of the effective fracture aperture. Normal 
stress changes are expressed relative to the acting normal equilibrium stress state which is equivalent 
to the sum of induced contact forces and the constant equilibrium normal stress state. Effective 
aperture fluctuations are induced by perturbations of the equilibrium state by means of fluid pressure 
changes or vice versa. Right: Consideration of the fracture’s microstructure in terms of a contact 
factor, respectively controlling the difference between initial hydraulic and mechanical opening 


High resolution XRCT characterization 


Fig. 4.3 Uniaxial setup for the determination of fracture stiffness E™ under harmonic excita- 
tion at different amplitudes in combination with fracture roughness characterization through high- 
resolution XRCT 


monolithic approaches (Schmidt and Steeb 2019). Staggered numerical schemes 
treat the fracture and the rock bulk domain individually choosing solvers for both 
parts independently. Thus, highly efficient solvers for poro-elastic equations and 
for pressure diffusion equations could be adopted. This is especially interesting if 
large problems need to be solved, e.g. on HPC platforms (Schmidt et al. 2022b). 
Communication between the fracture domain BF” and the the rock domain B?? 
is conducted via boundary conditions, e.g. through special software packages like 
preCICE. Constraints have to be introduced in order to guarantee numerical stability. 
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Solid Domain 


HD Domain 
(a) BFr (b) 


Fig. 4.4 Comparison of the weak coupling approach using non-conformal meshes (a) and a strong 
coupling approach with implemented interface elements (interface elements, and auxiliary nodes 
are explicitly shown for presentation purposes only; in the final discretization, the fracture surfaces 
align with d = 0) (b). Figure from Schmidt and Steeb (2019) 


In contrast, monolithic coupling approaches build a single algebraic set of equations 
which is solved for both domains simultaneously. 

Monolithic coupling of the fracture and the poro-elastic domain requires a math- 
ematical discussion of the balance equations introduced to govern hydro-mechanical 
flow processes within the fracture, cf. Schmidt (2022), Schmidt and Steeb (2019). 


4.2 Lower-Dimensional Representation of Frictional 
Interfaces 


For constructing the weak form of the quasi-static equilibrium conditions 
divo +og=0 (4.4) 
in the presence of displacement jumps we introduce an enriched test function v 


consisting of a continuous (standard) part v, and a Heaviside enrichment H (x)ar in 
the form 


v = V, + A(x)ar (4.5) 
with 
—0.5 Yxe Q 
H (x) = 4.6 
= los veeat eg 


Therefore, (vt — v~ )xr = [v]r = ar represents the jump in the test function across 
the embedded discontinuity surface F. Superscripts + and — denote the opposite 
sides of the discontinuity surface. Then, we find the gradient of the test function as 


grad v = grad v: + 6p(x)[v]r 8 n7 (4.7) 
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where the Dirac function is defined and normalized as 


Vv T 
op (X) = es with [ro dQ = 1 (4.8) 
x VxeT 
Q 
To construct the weak form of Eq. (4.4) 
0= |v [dive + og] dQ (4.9) 


Q 


we perform partial integration of the stress divergence term by employing 
divo -v = div (øv) —o: grad ve — 0 : p(x) v]r 8 n7 (4.10) 


which results in 


I [o : grad v — ob: v] d2+ [on -ivirar = ft- var (4.11) 
Q\r T IQ 


Withon” = t” and tt = tr = -t” we obtain 


f [o : grad v — ob: v] a2- f te-tvirar = [ t-var (4.12) 


Ar T IQ 


Note that o and tr are the total stresses in aHM formulation of a fluid-saturated 
porous medium and thus couple into the following flow equations by means of the 
effective stress principle: 


o =0'— appl (4.13) 
tr = tp — a, pnr (4.14) 


Note also that the present work assumes continuity of the pressure between matrix 
and fracture (p" = pf = p). 

In the above o’ is the effective stress tensor in the matrix, œg and af are coefficients 
weighting the pore-pressure contribution in the matrix and fracture to be defined 
suitably, p is the fluid pressure, I is the identity tensor, © = porr + (1 — ®)osr is 
the bulk density of the porous medium comprising the densities of the liquid and the 
solid, OLR and osr respectively, of porosity & and g is the gravitational acceleration. 
In the fracture with local unit normal nr, tp represents the effective stress vector. 

The effective quantities have to be supplied by constitutive equations. For the 
matrix, they are here given in incremental form 


do’ =C™: de®, (4.15) 
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where C™ is the material stiffness tensor, € = 5 ( gradu + grad Tu) is the linear 
strain tensor, €° its elastic part and u is the solid displacement vector. 

For the fracture effective stress, a connection is made constitutively to the dis- 
placement jump. Similar to the test function we specify u = u. + H(x)w(xr) such 
that (ut — u”),. = [u]r = w. Then, in analogy to the incremental relationship for 
the matrix effective stress tensor, the fracture effective traction vector follows incre- 
mentally from elastic increments in the fracture face relative displacement vector 


dti = K'dw° (4.16) 


where dw“ is the elastic part of the fracture relative displacement increment and 
Kf is the second-order stiffness tensor of the fracture comprising normal and shear 
components. 

For the following examples, the interface failure behaviour is of particular interest. 
Here, a Coulomb-type frictional interface was chosen 


dt’ = Kallu] = [Kar 8 nr) + Ks(1— nr Q nr)] dw* (4.17) 


with a split of the displacement jump analogous to the strain split used in elasto- 
plasticity w = w° + wP. 

Based on the work by Coulomb in 1773 (Davis and Selvadurai 2002), the frictional 
failure criterion can be written as 


te =c—o, tang (4.18) 


and describes a linear relationship between the limiting shear stress t; in a plane and 
the normal effective stress o, (tension positive) acting on that plane based on the 
cohesion c and friction angle @ of the fracture. The plane is in our case given by the 
fracture plane and a yield function f, describing the failure surface can be defined 
as follows 


fy =t—c+o, tang. (4.19) 


where T is the acting shear stress in the fracture plane. For all states of stresses, which 
are within this envelope, i.e. f, < 0, no failure is predicted. This region grows as 
the fracture normal stress increases. Once fy = 0 is reached the fracture deformes 
plastically according to the flow rule 


Ogy 
ot’ 


dw? = da (4.20) 


where the plastic potential is formulated to allow for non-associated flow: 


gy =T +0, tan y. (4.21) 
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Fig. 4.5 Flow through discrete fracture network, where fractures are represented as 2D elements 
in a 3D setting (Lee 2022). Barsch & Nagel. Technische Universtität Bergakademie Freiberg. CC 
BY SA 


where w is the dilatancy angle. Furthermore, the Kuhn-Tucker complementary con- 
ditions hold 
frA=0 AzO fp<0 (4.22) 


From the implicit integration of the material model the consistent stiffness matrix is 
passed back to the finite element assembly routines. 


4.3 Verification and Application Examples 


The lower-dimensional interface approach can be used to simulate flow in fracture 
networks. An example from OGS simulations is shown in Fig. 4.5. 

In this chapter, however, examples for mechanical and hydromechanical effects 
will be described briefly. 


4.3.1 Coulomb Model Verification 


To test the functionality of the Coulomb implementation, a compressive stress is 
created by compressing an element assembly separated by an interface vertically 
by a tenth of the initial fracture aperture, i.e. by 1 um, in a displacement-controlled 
manner. The top element is then sheared horizontally into the positive x direction by 
10 um, followed by shear unloading and reloading into the negative x direction by 
10 um, while the vertical displacement is held constant. The latter boundary condition 
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together with apositive dilatancy angle has the effect ofincreasing the fracture normal 
stress (accumulation of plastic normal fracture opening). Note that during this test 
the bottom boundary is fixed by a zero-displacement boundary condition whereas 
the boundaries at the left and right are freely movable. 

Shearing causes a linear deformation in the horizontal direction until the applied 
shear exceeds the yield strength of the material. After this point, the shear stress 
increases at a lower rate until u, |,—>m reaches its maximum at maximum displace- 
ment, which illustrates the apparent hardening during the elasto-plastic deformation 
process under increasing confinement associated with the dilatant material response. 
Note that the material itself is not hardening but ideally plastic. Afterwards, the 
unloading /reloading due to shear reduces the shear stress in the same linear fashion 
observed in the initial stage of shearing (elastic unloading). After exceeding the yield 
strength, the material follows the dilatancy-induced apparent hardening again until 
maximum load reversal. The transition from linear deformation to material hardening 
after exceeding the yield strength clearly illustrates the effect of the Coulomb model 
on the overall deformation process. Also, the correct loading-unloading behaviour 
could be shown. 


4.3.2 Soil-Structure Interfaces 


Soil structure interactions play an important role in geotechnics. They include earth 
pressure development behind rough retaining walls, foundation tractions, interactions 
between tunnel liners and rock masses, driven piles, and many others. 

Here, a simple laboratory test is chosen for illustrative purposes. Oedometer tests 
are a fundamental means of material characterization in geotechnical engineering. If 
the oedometer ring is compliant enough, radial stresses can be inferred in addition 
to the axial load in what is often called a soft oedometer. This example illustrates 
how the LIE model can be used to allow for relative displacement between the soil 
sample and the deformable walls of the soft oedometer. 

The soft oedometer setup consists of two domains. The larger one represents the 
soil sample and the smaller one the compliant metal ring of the oedometer. The 
top boundary is assigned a constant normal traction representing the applied load. 
The bottom of the oedometer is fixed in the vertical direction and free to move in 
horizontal direction. The top and the right side of the metal ring are defined as traction 
free. 

Due the Coulomb criterion, the soil can slide along the metal ring of the oedometer 
with a constant friction angle. Here, however, a frictionless interface was modelled. 
As illustration of the displacement discontinuity is the objective, the geomaterial is 
represented crudely oversimplified by a linear elastic model. The material properties 
of the benchmark are listed in the following table (Table 4.1): 

Under the action of the top load, the soil gets compacted vertically (Fig. 4.6). 
The metal ring on the other hand experiences nearly no vertical displacements. The 
displacement jump in the solution is clearly visible by the transition from the smooth 
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Table 4.1 Material properties 


Young’s modulus of soil Eı 150 MPa 
Young’s modulus of metal ring E2 210 GPa 
Poisson’s ratio of soil vı 0.2 

Poisson’s ratio of metal ring v2 0.3 

normal stiffness of the fissure Kn 1x10!5 Pa/m 
shear stiffness of the fissure Ks 0 Pa/m 

Top load p 1 MPa 
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Fig.4.6 Vertical displacements in the soft oedometer (OGS simulation). The use of frictional lower- 
dimensional elements allows for the establishment of a vertical displacement jump between the soil 
(smooth color gradient, left) and the oedometer ring (red zone with zero vertical displacements, 
right) in the solution on a single domain 


color gradient to the red outer part of the domain. This jump is inherent in the finite 
elements pace due to the Heaviside enrichment described above. 

For the chosen parameters and numerical settings, the radial displacements caus- 
ing the expansion of the ring used for radial stress measurements reach about 5.5 um 
at the circumference and cause radial stresses of around 0.229 MPa. This value is 
slightly lower than the 0.25 MPa that would be expected for a rigid ring. Likewise, 
the axial compaction of the sample itself is somewhat larger than in the case of a 
rigid ring. The exact values can be confirmed analytically by solving the linear elastic 
equilibrium problem for the sample subject to the top load boundary condition in 
conjunction with the pressure vessel equation, where both domains are connected 
via radial stress and displacement continuity. 


4.3.3 Fault Slip Experiment (Mont Terri) 


Simulations related to the FS experiment in Mont Terri (Guglielmi et al. 2017) 
are described elsewhere (Rutqvist et al. 2020). Here, we show more recent results 
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Fig. 4.7 3D model geometry 
for calculation of injection 
induced fault activation 


20 m 


using the LIE implementation described above. This example also links methods and 
application settings of WPs 2 and 3. 

To simulate the fault slip displacement induced by fluid injection-induced over- 
pressures, we define a 3D model domain (Fig. 4.7), which broadly represents the 
minor fault geometry embedded in the rock matrix (Rutqvist et al. 2020). The model 
domain has edge lengths of 20m and contains the fault dipping at 65°. Symmetry in 
the x direction is assumed and used to reduce the computational effort, therefore the 
side-length in the x direction reduces to 10m. 

The host rock matrix is modelled using linear elasticity assuming an isotropic 
material behavior with a bulk modulus of 5.9GPa and a shear modulus of 2.3 GPa, 
based on the average values derived from experiments on Opalinus Clay at Mont 
Terri. For all outer boundaries of the matrix we assume roller boundary conditions. 
Based on a simplified experimental data set, the in-situ stress field is characterized 
by oz; = —7 MPa ‚o,, = —6 MPa o,, = —3.3 MPa, where negative values indicate 
compression. Hydraulically, the host rock is considered impermeable because of the 
very low permeability of Opalinus Clay and the relatively short time-frame of these 
experiments. Fracture flow with various permeability models (constant, cubic law, 
cubic law only after shear slip, etc.) is applied while the fracture is assigned Coulomb 
behaviour. In this hydraulic-mechanical approach, a sudden increase in hydraulic 
aperture by a pre-set factor is modelled once shear failure occurs (Rutqvist et al. 
2020). After fault activation, further aperture changes due to hydraulic-mechanical 
coupling are calculated by the cubic law. To parameterize the fluid phase and the 
porous medium we use standard values, see Table 4.2. 

The initial fluid pressure estimated from site measurements was set to 0.5 MPa 
while the effective fracture traction is calculated from the stress field in the matrix 
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Table 4.2 Material properties for the modelling scenarios of the fault-slip experiment, scenario 1 


Material Parameter Value 
Scenario Nr. 
Fault Normal stiffness Kn [GPa/m] | 20 
Shear stiffness K, [GPa/m] 20 
Static friction angle [°] 22 
Tensile strength [MPa] 0 
Permeability model Cubic law after slip 
Host rock matrix Bulk modulus K [GPa] 5.9 
Shear modulus G [GPa] 2.3 
Bulk density @ [kg/m?] 2450 
Permeability [m7] 0 
Fluid Density o [kg/m?] 1000 
Compressibility [Pa~!] 4.4e-9 
Dynamic viscosity u [Pas ] 1.0e-3 
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Fig. 4.8 Injection volume and fracture aperture change (numerical and experimental) using fully- 
coupled LIE models driven by step-wise injection pressure increases 


by appropriate tensor rotation operations. The pressure and stress fields are assumed 
to be uniform and we neglect gravity in this case. 

An initial comparison between the results of the numerical calculation and the 
laboratory data still shows an offset in the results (Fig. 4.8). Currently, parameter 
studies on the initial stress field, the hydraulic and mechanical fault properties, the 
activation mechanism and the injection control are undertaken to improve the results. 
Other modelling approaches were published in an overview paper (Rutgvist et al. 
2020). Another modelling approach implemented in OGS was published in Urpi 
et al. (2020). These references contain details on the hydraulic models and the inter- 
pretation of the experimental data. A publication based on the approach presented 
here is currently under preparation. 
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Fig. 4.9 Area of fault affected by shear slip 


The area affected by shear slip, in which the fault transmissibility is drastically 
changed due to fault activation, is highlighted in Fig. 4.9. The area is centered around 
the injection point and is roughly circular in shape. This is inline with other modelling 
studies (Rutgvist et al. 2020). 


4.4 Concluding Remarks 


Lower-dimensional continua are mechanically and numerically intricate ways of 
representing interfaces of various kinds, including contact zones, fractures, faults, 
etc. Conceptually, the co-dimensional entities are generated by integration in their 
normal direction while making specific assumptions on the distribution of certain 
field quantities in this direction. Numerically, they can be represented by enrichment 
of finite element spaces and interface elements, for example. This kind of represen- 
tation comes with both advantages and disadvantages compared to other methods. 
While this chapter showed examples on hydraulic, mechanical and hydromechani- 
cal behaviour of elastic and frictional interfaces, further aspects related to fracture 
mechanics and hydraulic fracturing are discussed in Yoshioka et al. (2019), Mollaali 
et al. (2022). 


4 Pathways and Interfaces Under Stress Redistribution 75 


References 


Bandis, S., Lumsden, A., & Barton, N. (1983). Fundamentals of rock joint deformation. Interna- 
tional Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 20, 249-268. 

Biot, M. A. (1941). General theory of three-dimensional consolidation. Journal of Applied Physics, 
12(2), 155-164. 

Davis, R. O., & Selvadurai, A. P. S. (2002). Plasticity and geomechanics. Cambridge: Cambridge 
University Press. 

Gens, A., Carol, I., & Alonso, E. (1990). A constitutive model for rock joints formulation and 
numerical implementation. Computers and Geotechnics, 9, 3-20. 

Guglielmi, Y., Birkholzer, J., Rutqvist, J., Jeanne, P., & Nussbaum, C. (2017). Can fault leakage 
occur before or without reactivation? Results from an in situ fault reactivation experiment at Mont 
Terri. Energy Procedia, 114, July 2017, 3167-3174. 

Lee, H. (2022). Solving the rock-hard problem of nuclear waste disposal burial is best. Ars Technica, 
1-11. 

Mollaali, M., Kolditz, O., Mengsu, H., Park, C.-H., Park, J.-W., McDermott, C., Chittenden, N., 
Bond, A., Yoon, J., Zhou, J., Pan, P.-Z., Liu, H., Hou, W., Lei, H., Zhang, L., Nagel, T., Barsch, 
M., Wang, W., Nguyen, S., Kwon, S., Lee, C., & Yoshioka, K. (2022). Comparative verification 
of hydro-mechanical fracture behavior: Task g of international research project DECOVALEX 
2023. International Journal of Rock Mechanics and Mining Sciences (under review). 

Renner, J. (2020). STIMTEC: A mine-scale hydraulic stimulation experiment of anisotropic meta- 
morphic rock with evaluation by mine-back drilling. ARMA Newsletter. 

Rutqvist, J., Graupner, B., Guglielmi, Y., Kim, T., Maßmann, J., Nguyen, T. S., et al. (2020). An 
international model comparison study of controlled fault activation experiments in argillaceous 
claystone at the Mont Terri laboratory. International Journal of Rock Mechanics and Mining 
Sciences, 136, 104505. 

Schmidt, P. (2022). Hydro-mechanical coupling of flow in deformable high-aspect ratio fractures. 
Ph.D. thesis, University of Stuttgart. 

Schmidt, P., & Steeb, H. (2019). Numerical aspects of hydro-mechanical coupling of fluid-filled 
fractures using hybrid-dimensional element formulations and non-conformal meshes. GEM - 
International Journal on Geomathematics, 10(1), 1-36. 

Schmidt, P., Steeb, H., & Renner, J. (2021). Investigations into the opening of fractures during 
hydraulic testing using a hybrid-dimensional flow formulation. Environmental Earth Sciences, 
80, 497. 

Schmidt, P., Dutler, N., & Steeb, H. (2022a). Importance of fracture deformation throughout 
hydraulic testing under in situ conditions. Geophysical Journal International, 228, 493-509. 
Schmidt, P., Jaust, A., Steeb, H., & Schulte, M. (2022b). Simulation of flow in deformable fractures 
using a quasi-newton based partitioned coupling approach. Computational Geosciences, 26, 381- 

400. 

Steeb, H., & Renner, J. (2019). Mechanics of poro-elastic media: A review with emphasis on 
foundational state variables. Transport in Porous Media, 130(2), 437-461. 

Urpi, L., Graupner, B., Wang, W., Nagel, T., & Rinaldi, A. P. (2020). Hydro-mechanical fault 
reactivation modeling based on elasto-plasticity with embedded weakness planes. Journal of 
Rock Mechanics and Geotechnical Engineering, 12(4), 877-885. 

Vinci, C., Renner, J., & Steeb, H. (2014). A hybrid-dimensional approach for an efficient numerical 
modeling of the hydro-mechanics of fracture. Water Resources Research, 50, 1616-1635. 

Watanabe, N., Wang, W., Taron, J., Görke, U. J., & Kolditz, O. (2012). Lower-dimensional interface 
elements with local enrichment: Application to coupled hydro-mechanical problems in discretely 
fractured porous media. International Journal for Numerical Methods in Engineering, 90(8), 
1010-1034. 


76 M. Barsch et al. 


Yoshioka, K., Parisio, F., Naumov, D., Lu, R., Kolditz, O., & Nagel, T. (2019). Comparative verifi- 
cation of discrete and smeared numerical approaches for the simulation of hydraulic fracturing. 
GEM - International Journal on Geomathematics, 10(1), 13. 


Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 
International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, 
adaptation, distribution and reproduction in any medium or format, as long as you give appropriate 
credit to the original author(s) and the source, provide a link to the Creative Commons license and 
indicate if changes were made. 

The images or other third party material in this chapter are included in the chapter’s Creative 
Commons license, unless indicated otherwise in a credit line to the material. If material is not 
included in the chapter’s Creative Commons license and your intended use is not permitted by 
statutory regulation or exceeds the permitted use, you will need to obtain permission directly from 
the copyright holder. 


Chapter 5 A) 
Virtual Reality and Computational get 
Efficiency 


Karsten Rink, Nico Graebling, Lars Bilke, Jörg Buchwald, Thomas Fischer, 
Christoph Lehmann, Tobias Meisel, Dmitri Naumov, Wenqing Wang, 
Keita Yoshioka, and Olaf Kolditz 


In this chapter we briefly describe information methods and technologies support- 
ing geotechnical systems analyses, i.e. using virtual reality methods for data and 
model integration (Sect. 5.1) and improving computational efficiency by using high- 
performance-computing techniques (Sect. 5.2). 

The generalised approach of the developed visualisation methodology makes it 
possible to integrate a variety of data formats into 3D scenes and to create studies 
for a wide range of application fields! For example, a prototype for an experiment 
information system was developed for the rock laboratory in Mont Terri operated 
by the Federal Office of Topography surveys Switzerland (swisstopo) (Bossart et al. 
2017; Jaeggi et al. 2017). This prototype combines a representation of the complex 
geometry of the Swiss Jura mountains with tunnel system and boreholes of the rock 
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laboratory with a series of results of coupled numerical simulations for the planning 
and validation of long-time experiments, some of which are designed to run over 
several decades. 

Due to the resolution, runtime and complexity of the simulations, the previously 
described methods for data reduction are essential for the visualisation of the results 
(Fig. 5.3). Figure 5.1 depicts amap ofthe URL Mont Terri including local labels for 
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the experiments which have been integrated into the Mont Terri visualization study. 
The development was carried out in close cooperation with swisstopo as well as the 
expert scientists who are in charge of the experiments. 

The visualisation of geotechnical processes requires three-dimensional represen- 
tation options for coupled field problems as well as validation in the context of 
corresponding measurement results, therefore, in the following we discuss aspects 
of data and model integration into a VR framework. 


5.1 Visual Data and Model Integration 


A prototype of a Digital Twin for URL Mont Terri has been designed and imple- 
mented (Graebling et al. 2022). In Sect. 2.6 the importance for data visualization and 
placing models into the real geological context has been already elaborated. In order 
to set-up the Digital Twin a particular Task VR (virtual reality) has been established 
in the portfolio of the Mont Terri project which was supported by the consortium. 
The VR task is dedicated to both data and model integration in order to combine the 
two sources of information for process characterization in the subsurface. The VR 
task has been evolved as a starting point of a new kind of data management for the 
future-a visually mediated data base including predictive capacities by integrated 
process models. 

From the technical side, several software tools are important for the development 
ofthe VR study. The DataExplorer as part of OpenGeoSys-Workflows is used for data 
integration (Rink et al. 2022). The Unity framework is being adopted for interactive 
visualization purposes and data base access interfaces have been implemented. 


5.1.1 Data Integration 


An important prerequisite for data integration was linking two data bases, (i) the 
borehole information system BIS (developed a hosted by Simultec) and (ii) the 
sensor data base Geoscope (developed and hosted by Sixense Soldata). The support 
from the two companies in helping for the data base links is greatly acknowledged. 

An illustration of data integration for URL Mont Terri is depicted in the figure 
collection. Figure 5.2a shows the entire tunnel system of URL Mont Terri includ- 
ing the positions of all active experiments (green spheres). The acronyms stand 
for the respective experiments. The indicated positions of the experiments can be 
clicked in the Unity framework and further information for connected boreholes and 
experiments become available. By clicking a particular borehole related geometri- 
cal information is displayed in a borehole menu (Fig 5.2b). Measurements such as 
displacements due to rock convergence can be mapped on the related boreholes of 
an experiment, therefore, structural and process data are being combined (Fig. 5.2c). 
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Fig. 5.2 Data integration for URL Mont Terri: Prototype CD-A experiment 


Structural information such as fracture planes from surface scans can be added into 
the scene (Fig. 5.2d). Other documentation such as related reports, publication links 
or photos from cores can be displayed in the context of the selected experiment or 
borehole (Fig. 5.2e). Moreover, Fig. 5.2f shows sensor data for the temporal evolution 
of displacements measured from different sensor positions. Therefore, the experi- 
mental information systems allows to show and combine structural and process data 
distributed in time and over space. 
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5.1.2 Model Integration 


Model integration of three examples of Mont Terri experiments are presented: 


e CD-A: Influence of Humidity on the Cyclic and Long-Term Deformation Behavior 
ofthe Opalinus Clay: The CD-A experiment focuses on the evaluation of the highly 
coupled hydraulic-mechanical (HM) processes, aiming on the identification of 
the significant mechanical effects in the vicinity of an excavation. To gain an 
increased process understanding of the excavation induced mechanical processes, 
a comparison of coupled hydraulic-mechanical effects in the CD-A twin niches— 
one containing coupled hydraulic-mechanical effects due to desaturation and the 
other one characterized by avoiding desaturation effects as much as possible, is 
carried out (Ziefle et al. 2022). 

FE: The “Full-Scale Emplacement” experiment mimics, as realistically as possible, 
the construction, waste emplacement, backfilling and early post-closure evolution 
of high-level waste disposal in Opalinus clay by a multiple heater test. It aims at 
the investigation of repository-induced thermo-hydro-mechanical (THM) coupled 
effects on the host rock at this scale and the validation of existing coupled THM 
models (Miiller et al. 2017; Papafotiou et al. 2019) 

FS-B: The “Fault Slip” experiment studies the aseismic-to-seismic activation of a 
fault zone in a clay/shale formation of Opalinus clay and particularly the condi- 
tions for slip activation and stability of faults. Furthermore, implications of fault 
permeability and the coupled effects between fault slip, pore pressure are investi- 
gated (Nussbaum et al. 2017; Guglielmi et al. 2017, 2020; Shiu et al. 2021; Cappa 
et al. 2022). 


Simulation results for the three experiments using OpenGeoSys are integrated into 
the VR study showing the process evolution in the real context on the URL Mont 
Terri.” In case of the CD-A experiment, the numerical simulation has been used for 
the experimental design (calculation of optimal distance between the two niches) 
(Fig. 5.3). 


5.2 Computational Efficiency 


In order to achieve sufficient computational efficiency for the simulation of strongly 
coupled and often non-linear multi-field problems such as for thermo-hydro- 
mechanical (THM) processes, massive parallelization technologies have to be applied. 

OpenGeoSys (OGS) is an open source software developed at the UFZ for the sim- 
ulation of coupled thermal-hydrological-mechanical-chemical (THMC) multi-field 
processes in porous media. Such computations place high demands on memory and 
computing power. The aim is to model real-world applications by increasing spatial 


? The Virtual Reality study for the Underground Research Laboratory Mont Terri (VR Task) has 
been supported by the Mont Terri consortium which is greatly acknowledged. 
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Fig. 5.3 Model integration for selected experiments for URL Mont Terri 


resolutions and to make models more realistic, for example by integrating non-linear 
material laws into the modelling. With a high spatial resolution, important geological 
structures that are relevant for the physical processes in the domain can be incorpo- 
rated into the simulation model. This allows, for example, detailed simulations in 
large-scale areas as they are necessary for predictions about the mechanical integrity 
of various host rock types. The implementation required an efficient use of the avail- 
able HPC resources. Within the framework of the project, scalability and parallel 
efficiency of individual OGS components (e.g. I/O, various local assemblers, mono- 
lithic approach vs. staggered scheme, linear equation system solvers) have been 
analysed. For the HPC optimisation, tracing/profiling tools were used to identify 
time-critical code sections and to find improvement potentials. Possible re-design 
steps for OGS were then derived from the analysis results and implemented. The 
optimised OGS code was successfully demonstrated in various application areas, 
such as hydrology and environmental geotechnics. 

One of these time-critical code section was identified as the simulation output file 
operations. Before we used parallel file I/O based on the Visualization Toolkit (VTK) 
but its file formats were not designed with massive parallel clusters and file systems 
in mind and were a serious bottleneck in some application scenarios. We therefore 
implemented simulation result output based on HDF5 which decreased the time for 
writing simulation output files by more than an order of magnitude. Post-processing 
can still be done with common visualization tools like Para View. 

The variational phase-field model in this study has been implemented in parallel 
in OpenGeoSys. Unlike numerical methods that treat fracture explicitly with element 
edge or node splitting, the variational phase-field model treats fractures through a 
phase-field variable, which is a primary variable to be solved. Therefore, its parallel 
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Fig. 5.4 THM simulations for a complete deep geological repository, shown is the temperature 
distribution 100 years after emplacement (Wang et al. 2021) (left), THM simulations for the full- 
emplacement (FE) experiment in the Mont Terri underground laboratory, shown are temperature, 
saturation and vertical stress distributions (Raith et al., 2020) (right) 


implementation does not require any special consideration. If one has an implemen- 
tation of linear elasticity, all that is required is to add an extra equation system for 
phase-field evolution and iterate the linear elastic system and the phase-field evo- 
lution system. All the variational phase-field simulations in this study were run in 
parallel and the linear scaling has been confirmed up to 100’s of CPUs. 

The application examples depicted in Fig. 5.4 present complex thermo-hydro- 
mechanical (THM) simulations for geotechnical problems such as the design of deep 
geological repositories for the storage of energy waste (e.g. heat-generating radioac- 
tive substances) in the geological subsurface. The HPC methodology is particularly 
demanding here because, on the one hand, it deals with non-linear coupled multi-field 
problems that place high demands on the parallel solution methods. Secondly, they 
involve complex geometries (mapping of geological and geotechnical structures), 
which place high demands on the domain discretisations and corresponding domain 
decompositions. 


5.3 Software Engineering and Benchmark Workflows 


5.3.1 Containers 


Linux containers are a technology for encapsulating the runtime environment of 
software and executing it on any Linux system. Common technologies here are 
Docker and Singularity. Singularity is ideal for HPC environments because it enables 
seamless integration into parallelised environments (e.g. with MPI) and also uses a 
better execution concept than Docker (Singularity containers can be executed on 
the system with limited user rights whereas Docker containers require administrator 
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Fig. 5.5 Schematic workflow for container generation 


rights). Singularity is thus widely available on HPC systems, for example on the 
systems of the TU Dresden (“Taurus”), the Forschungszentrum Jülich (“Juwels”) 
and the Helmholtz Centre for Environmental Research (“Eve”). 

In order to provide various software configurations of OpenGeoSys (e.g. different 
OGS versions, integration of third-party libraries to extend the range of functions) 
through Linux containers, a container creation framework “OpenGeoSys Container 
Maker” was implemented, which provides a user-friendly parameterisation of the 
container to be created. The framework is available to users via a web interface 
on the OpenGeoSys GitLab server (used for software version management, project 
management and software quality assurance). The container is created according to 
user specifications and made available as a download (Fig. 5.5). 

The containers include the OpenGeoSys software in binary form as well as other 
tools helpful for the preparation and execution of simulations, such as computational 
grid generators or Python scripting interfaces. Singularity containers are regular files 
and can be transferred onto the HPC system easily after creation via the web interface 
mentioned above. Parallelised simulations can now be carried out with the help of 
the container. We used these container-based OGS workflows successfully during 
the project. 


5.3.2 Benchmarks and Jupyter Notebooks 


A new development within the project is the provision of container applications for 
Jupyter Notebooks.” In addition to the OGS core and external modules / libraries (e.g. 
MFront, PHREEQC, PETSc), these containers also contain the Jupyter Notebook 
server application and a number of Python packages which can be added as needed 
(Fig. 5.6). After starting the container, the Jupyter Notebook can be accessed as a 
browser application and OGS can be executed using notebooks. Jupyter Notebooks 
also form a new basis for OGS benchmark presentation and integration. New test 
cases are formulated and explained in Python-based Jupyter notebooks which can 
intermix script logic with explanatory text and images. Moreover, the large variety 
of existing Python tools can be used for pre- and postprocessing of OGS simulation 


3 https://www.opengeosys.org/docs/userguide/basics/jupyter-notebooks/. 
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results. Figure 5.8 shows the OGS benchmark gallery page* which is organized 
according to the THMC process coupling hierarchy (Fig. 5.7). 


5.3.3 Benchmark Workflows 


The OpenGeoSys benchmark gallery (Fig. 5.8) is organized according to the 
THM/RTP process hierarchy, thermo-hydro-mechanical and reactive transport pro- 
cesses. A specific process class is represented by a tile showing a simulation result of 
the related process class. After clicking a tile, the available benchmarks of a process 
class comes visible and can be selected. Typically, an OGS benchmark starts with a 
short description of the problem and showing most important results for the bench- 
mark test. All benchmarks are linked with the OGS project file (prj-file), therefore, 
the benchmark settings are directly available through the gallery. All benchmarks 
are part of the OGS quality assurance workflow which is continuously running all 
tests (benchmarks and so-called unit-tests for basic functionalities) after any code 
changes, automatically. For new benchmarks Jupyter notebooks are available for 
user convenience and user-specific pre- and postprocessing operations. New process 


4 https://www.opengeosys.org/docs/benchmarks/. 
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Fig. 5.8 OGS Benchmark Gallery organized by process classes. Benchmarks of a specific process 
class are behind the tiles 


classes and/or those with new benchmarks are highlighted as featured processes on 
top of the benchmark gallery. 

Big Data applications cover a wide range of requirements for the analysis of 
complex data. This usually includes not only the analysis itself, but also the nec- 
essary pre-processing steps for data integration and preparation as well as methods 
for user interaction and evaluation. Therefore, it is necessary to create and execute 
complex workflows even for seemingly simple analyses. Modern HPC architectures 
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are an ideal basis for providing customised and therefore high-performance working 
environments for different application requirements. 

Large amounts of data not only reach the limits of computability, but above all 
of comprehensibility. Simulation processes are usually divided into three parts: in 
pre-processing, the simulation model is defined and input data for the simulation 
is created. Typically, the simulation as the second part generates large amounts of 
result data, which must then be processed and analysed in post-processing. Post- 
processing is usually based on a visualisation of the result data and takes place 
either on front-end computers of the simulation cluster or on the user’s desktop 
computer. Only in the latter case a transfer of the simulation data is necessary. The 
final analysed and processed data as a result of post-processing are far less extensive 
than the simulation results. Performing post-processing directly during the simulation 
(in-situ visualisation) makes it possible to bypass the communication/data transfer 
bottleneck and transfer only the already analysed data to the user’s desktop computer. 

The development, numerical realisation, validation and use of integrative mod- 
elling tools for the simulation of problem-specific, coupled processes as well as 
the efficient data presentation in the Big Data context make essential methodolog- 
ical contributions to the systematic analysis and prediction of processes in various 
areas of environmental science, especially with temporally and spatially large model 
dimensions. In the overarching conceptual and methodological approach to model 
and software development, the project results achieved by the grant recipient UFZ 
play a leading scientific role in the efficient and sustainable planning and management 
of the environmental science systems under investigation. 

Work on the further development of the OpenGeoSys software platform was nec- 
essary in order to be able to reliably and efficiently carry out simulations of complex 
coupled processes of environmental science applications (hydrology, geotechnics) 
on large scales using HPC systems. This provided evidence that OpenGeoSys is suit- 
able for modelling environmental science processes at real sites. In this context, the 
research work on parallelising the software should be mentioned as necessary and 
appropriate. With regard to scientific 3D visualisation, the work carried out on the 
development of workflows, numerical algorithms and software modules for the syn- 
optic in-situ representation of heterogeneous model data from conceptually different 
sources has laid the foundations for the integrated consideration of structural infor- 
mation, results from measurement campaigns and results of numerical simulation, 
which were not available at the beginning of the project work.° 


5 OGS workflow development and High-Performance-Computing (HPC) capabilities have been 
supported by the ScaDS.AI and Pilot-Lab ESM initiatives which is greatly acknowledged. 
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Chapter 6 A) 
Synthesis and Outlook get 


Olaf Kolditz, Tuanny Cajuhi, Ralf-Michael Günther, Holger Steeb, 
Frank Wuttke, Keita Yoshioka, Norbert Grunwald, and Thomas Nagel 


The principal interest of the GeomInt project consists of the investigation of effects 
on barrier integrity of three host rock formations: clay, salt and crystalline. The 
project focuses on distinct physical processes that can influence barrier integrity 
in these rocks, particularly those related to swelling and shrinkage, pressure-driven 
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percolation and stress redistribution. In the first part of the project, methodologies 
were developed in both experimental and numerical areas (Kolditz et al. 2021b). 
Combined model-experiment studies (MEX) were used as the first synthesis tool. 
The problem classes (processes and rock types) were systematized and model and 
code comparisons were made for test examples (benchmarks) and laboratory exper- 
iments. A further result of the first project phase was also to systematically evaluate 
(advantages and disadvantages) the suitability of certain experiments and numerical 
methods for the individual or combined process classes (see above), which can influ- 
ence the barrier effect, and to make further studies available (Kolditz et al. 2021a). In 
the second phase of the project, the methodical tools were used in particular for the 
analysis of laboratory and field experiments. Model chains (Sect. 6.1) and formation- 
specific workflows (Sect. 6.3) were developed as further synthesis tools, which are 
briefly presented below. Moreover, the GeomInt project has benefited from the out- 
reach to various national and international activities such as the Mont Terri project, 
BenVaSim, DECOVALEX, and EURAD (Sect. 6.2). 


6.1 Mechanical Integrity of Clay Rocks 


6.1.1 Model Chains 


Barrier integrity can be affected by a number of processes. In the case of gas per- 
meation, relevant processes are highlighted in Fig. 6.1. Gases can arise e.g. from 
corrosion, degradation of organic compounds, setting of cement or pyrolysis. Gases 
initially spread out in dissolved form diffusely and advectively. At higher pressures, 
a gaseous phase can form, which propagates further as a multiphase flow. If the gas 
pressure continues to rise, mechanical damage to the rock with dilatancy and eventu- 
ally propagating fractures can occur. The EURAD project is currently investigating 
whether this damage in clay rocks is reversible when the fluid pressures decrease and 
the plastic rock material converges and the fractures close again (the so-called self- 
sealing effect). This process chain shown in Fig. 6.1 can be mapped using new calcu- 
lation methods (Grunwald et al. 2022). Dis- and continuum mechanical approaches, 
which were developed in GeomInt, are combined. The variational phase field (VPF) 
method is used for simulating fracture propagation processes (Yoshioka et al. 2022). 


6.1.2 Benchmarking 


For a reliable description of complex thermo-hydro-mechanical (THM) and reactive 
transport (RTP) processes, extensive tests of both the physical/chemical plausibility 
and the correct implementation in numerical simulators are essential. Figure 6.2 
shows an sketch of benchmarks for THM/RTP models. 
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Fig. 6.1 Model chain for simulating consecutive processes of gas propagation in clay rocks leading 
to damage barrier integrity functions (Marschall et al. 2005), extended by Grunwald 


Fig. 6.2 Sketch of 


=. 
benchmark arrangement of Tht \ S 
THM/RTP models T (BE 

| 4c | (= izae | 


An extensive suite of test examples was created for the development of the OGS- 
TH2M model (Grunwald et al. 2022). A hierarchical concept was developed and 
implemented that systematically checks all (reasonable) process couplings based on 
the individual processes (Fig. 6.3). The extensive collection of benchmarks (> 100 
tested test examples) from OGS is directly available via the portal and offers users 
an ideal introduction to THM/CB modelling. Some of the OGS benchmarks are 
already available as Jupyter notebooks (Sect.5.3.2) and can thus be integrated into 
other Python applications with the corresponding advantages of user-specific data 
analysis. 

The OGS-TH2M model, a new process class, developed by Grunwald et al. (2022). 
A comprehensive benchmarking has been conducted with the BenVaSim (Lux et al. 
2021; Pitz et al. 2022) and DECOVALEX projects (Sect. 6.2.1). 
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Fig. 6.3 Benchmarking hierarchy for TH2M processes (Grunwald et al. 2022) 


6.2 International Collaboration 


Fundamental developments from the GeomInt projects are widely used and deployed 
in various international collaboration activities in particular DECOVALEX, EURAD 
and the Mont Terri projects. 


6.2.1 DECOVALEX 


DECOVALEX (DEvelopment of COupled models and their VALidation against 
EXperiments) is an international research and model comparison project, initi- 
ated in 1992, for advancing the understanding and modeling of coupled thermo- 
hydro-mechanical-chemical (THMC) processes in geological and geotechnical sys- 
tems. DECOVALEX is an international community effort. OpenGeoSys is actively 
involved into all tasks of the current D2023 phase.' Task G is particularly dealing 
with the integrity of barrier rocks. This task is aiming to better understand reactivation 
of pre-existing discontinuities for brittle host rocks. In particular, the potential for 
existing features to undergo shear displacements and related changes in permeabil- 
ity as the result of coupled thermal, mechanical, hydraulic and chemical effects can 
all have significant impacts on repository safety functions (e.g., creating permeable 
pathways or, for very large displacements, mechanical damage of waste packages). 

Task G is organised in four steps with increasing complexity (Fig. 6.4). The 
Task route starts from fracture mechanics (Step 1: M process) then diverting into 
hydro-mechanical (HM) and thermo-mechanical (TM) processes and finally com- 
bined into thermo-hydro-mechanical (THM) processes. These steps are described in 
the following including the concept, experimental design and data which will form 


l https://decovalex.org/D-2023/overview.html. 
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Fig. 6.4 Task G structure En Step2 
with related steps ->j “ HM Stress-dependent 
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the foundation for benchmarking exercises and experimental analysis. The research 
results are presented in the interim DECOVALEX report and the essence in Mol- 
laali et al. (2022) which contains a comprehensive model and code comparison study 
including dis- and continuum mechanical approaches from all Task G members (nine 
international teams). Task G has been working on interpreting results from GeomInt 
laboratory experiments in Freiberg (Frühwirt et al. 2021) and Edinburgh (GREAT 
cell facility) (McDermott et al. 2018; Fraser-Harris et al. 2020). 


Step 4: 
Creativity part 
(TBD) 


S 


6.2.2 EURAD 


EURAD is a joint European project heading to safe radioactive waste management 
(RWM) through the development of a robust and sustained science, technology and 
knowledge management programme. It supports timely implementation of RWM 
activities and serves to foster mutual understanding and trust between Joint Pro- 
gramme participants.? GeomInt closely cooperated within EURAD particularly by 
contributing modeling know-how to work packages DONUT and GAS which are 
aiming at improving the understanding barrier integrity of clay formations in Europe 
suited as host rocks for deep geological repositories. GeomInt, therefore, actively 
helps in building European competences in clay systems understanding. Through 
EURAD national and international activities such as GeomInt, iCROSS, BenVaSim 
(Lux et al. 2021), DECOVALEX are bridged (Fig. 6.5). 


6.3 Workflows—Mont Terri Project 


One of GeomInt’s main synthesis services is the development and implementation 
of workflows for the analysis of geotechnical systems. This is to be shown as an 
example for the joint work in the Mont Terri rock laboratory. Only the main features 


? https://www.ejp-eurad.eu/. 
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eurad RTP 


GEOMINT 


Fig. 6.5 International collaboration scheme—process-based cooperation for building European 
competences 


of the workflow are shown here and the full description in Chap. 2 is referred to for 
details. (link Project geography) 

The experimental basis is based on BGR field experiments in the Mont Terri rock 
laboratory and subsequent laboratory experiments. For example, drillcore samples 
were taken from the sandy facies of the OPA and analyzed in the geotechnical lab- 
oratory of the University of Kiel. The anisotropic, fracture and shrinkage behaviors 
were examined (observation/monitoring) in detail. The CD-A field experiment at 
Mont Terri has played a central role as a bridge between laboratory and field scales 
(ref). Much of the data from the URL Mont Terri was visually integrated into an infor- 
mation system as described in (ref) (Data Integration) The laboratory experiments 
(e.g. on the fracture behavior of the clay samples) were simulated using different 
model approaches. The essential processes and parameters that are relevant for the 
barrier integrity could be identified. (Process Simulation/Data Analytics) On this 
basis, a first transfer of the models to the field scale could be carried out. Exemplary, 
we were able to model cracks due to excavation and due to desiccation using dis- 
crete and diffuse approaches. In addition, an experimental design was carried out 
in advance for the planning of the CD-A experiment: the optimal distance between 
the two niches of the CD-A experiment was determined by numerical modeling and 
then implemented accordingly. This is an example of the use of models for the plan- 
ning of complex and often expensive field experiments (Data Analytics/Knowledge 
Transfer). Finally, much of the experimental and modeled data and results have been 
integrated into the Visual Information System for Mont Terri and are available for 
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Fig. 6.6 Workflow concept applied to clay rock systems of the rock research laboratory in Mont 
Terri 


various planning, discussion and demonstration purposes (knowledge transfer) (Fig. 
6.6). 

In addition to the CD-A experiment, currently more experiments from the Mont 
Terri rock laboratory are integrated into the workflow concept, e.g. the FE (full 
scale emplacement heater) and FS (fault slip) experiments targeting at the “proof-of- 
concept” for comprehensive workflows for entire research labs. The laboratory tests 
carried out in this phase of the GeomInt project, e.g. fluid percolation, will serve as 
basis for the application of numerical methods to model in the latter in-situ tests. 
Future work will be dedicated to further professionalize the workflow concept and 
transfer the know-how and software framework to other research labs. 
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